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Abstract 

A formal derivation of linear hydrodynamics for a granular fluid is given. The linear response to 
small spatial perturbations of the homogeneous reference state is studied in detail using methods of 
non-equilibrium statistical mechanics. A transport matrix for macroscopic excitations in the fluid is 
defined in terms of the response functions. An expansion in the wavevector to second order allows 
identification of all phenomenological susceptibilities and transport coefficients through Navier- 
Stokes order, in terms of appropriate time correlation functions. The transport coefficients in 
this representation are the generalization to granular fluids of the familiar Helfand and Green- 
Kubo relations for normal fluids. The analysis applies to a wide range of collision rules. Important 
differences in both the analysis and results from those for normal fluids are identified and discussed. 
A scaling limit is described corresponding to the conditions under which idealized inelastic hard 
sphere models can apply. Further details and interpretation are provided in the paper following 
this one, by specialization to the case of smooth, inelastic hard spheres with constant coefficient of 
restitution. 
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I. INTRODUCTION 



Forty years ago, significant advances in the theory and simulation of simple atomic fluids 
were stimulated by the application of exact methods from non-equilibrium statistical me- 
chanics, namely linear response and the "time correlation function method" [1]. The results 
differed from earlier studies based on approximate kinetic theories in that they are formally 
exact and closely related to properties measured in experiments. Subsequently, a great deal 
has been learned through the study of appropriate time correlation functions by theory, 
simulation, and experiment [2]. In many respects, the more recent study of granular fluids is 
poised to exploit this body of work on normal fluids. Significant advances have been made 
in the past decade through the application of molecular dynamics simulation and kinetic 
theory. However, although the generalization of the formal structure for non-equilibrium 
statistical mechanics has been described [3, 4], relatively few applications outside of sim- 
ulation and kinetic theory [5-7] have been given. In particular, formally exact relations 
between properties of interest and appropriate time correlation functions appear restricted 
to the simplest cases of tagged particle motion [8-11] and liquid structure [12]. 

The objective here is to provide a general application of these formal methods to granular 
fluids. This is done by studying the response to small perturbations of a reference homo- 
geneous state, and using this to extract formally exact expressions for the hydrodynamic 
transport coefficients up through Navier-Stokes order. This is the analogue of the study 
of excitations about the equilibrium Gibbs state for normal fluids. In many experimental 
conditions of interest, for both normal and granular fluids, the system is not close to a 
global homogeneous state. Nevertheless, it is expected that the reference states studied here 
are relevant locally for more complex and realistic physical conditions [13]. For example, 
the transport coefficients such as viscosity and thermal conductivity obtained from linear 
response, are the same functions of density and temperature as those in the associated non- 
linear equations applicable under more general conditions. More specifically, the transport 
coefficients discussed here are the same as those studied to date for granular fluids using 
kinetic theory, but without the limitations or assumptions of those theories. 

The problem of linear response can be formulated at both the level of phenomenological 
hydrodynamics and statistical mechanics. In both cases, the response of the hydrodynamic 
fields y a (r,t) (local number density, granular temperature, and flow velocity), defined as 
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functions of position r and time t, to small initial spatial deviations 5y a (r,Q) from their 
values in a homogeneous reference state, is written in terms of a matrix of response functions 
C a/ 3 (r;t), 




These response functions can be calculated approximately using the linearized phenomeno- 
logical hydrodynamic equations, to obtain an explicit result as a function of the parameters 
(e.g., cooling rate, pressure, transport coefficients) in those equations. On the other hand, 
C a /3 (v, t) can be given a formally exact representation from non-equilibrium statistical me- 
chanics in terms of time correlation functions evaluated in the reference homogeneous state. 
If the hydrodynamic equations are valid on some length and time scale, then the exact and 
phenomenological representations must be the same in that context [14-17]. This provides 
a means for identifying the phenomenological parameters such as transport coefficients in 
terms of exact properties of the correlation functions, that is a precise link between the 
macroscopic properties of interest under experimental conditions and the underlying funda- 
mental microscopic laws. 

There are two parts to this prescription. The first is a demonstration that the statistical 
mechanical representation admits a limit with the same form as that from hydrodynamics, 
allowing identification of the hydrodynamic parameters. The second is a proof that the hy- 
drodynamics dominates all other possible excitations in this limit. The first part constitutes 
the usual derivation of Helfand [16] and Green-Kubo [17] time correlation function expres- 
sions for transport coefficients in a normal fluid, and the presentation here is essentially its 
extension to granular fluids. The second part is more difficult and remains incomplete even 
for normal fluids. The argument there is that the hydrodynamic fields have a dynamics that 
persists on the longest time scale, since they are the densities of conserved quantities and 
therefore have infinite relaxation times at infinite wavelengths. Hence it should be possible 
to wait for all other excitations to decay, leaving a space and time scale on which only hydro- 
dynamics remains. This is the limit in which the Green-Kubo expressions can be identified 
and hydrodynamics dominates. A similar assumption is made here for granular fluids, as 
discussed further below. 

The expressions obtained in this way for the phenomenological hydrodynamic parame- 
ters are still formal, in the sense that the time correlation function expressions are difficult 
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to evaluate. However, they depend only on the low frequency, long wavelength limit of 
the response functions and, therefore, do not entail the full complexity of evaluating the 
complete functions C a p (r; t). This approach gives directly exact expressions for the param- 
eters, without the intermediate practical assumptions and calculations needed to determine 
C a/ 3 (r; t) first. For example, it is entirely possible that the usual forms of kinetic theory, 
and the consequent derivation of hydrodynamics from them, might fail for granular fluids 
under some conditions. Still, hydrodynamics might apply independent of the kinetic theory 
basis and the correct expressions for the transport coefficients would be those obtained from 
linear response. The utility of such formal results for both simulation and theory, has been 
illustrated recently for granular fluids in the cases of mobility [10] and impurity diffusion 
[11]. 

Although the above prescription for application of linear response to the hydrodynamics 
of a granular fluid is simple to state in general, its implementation in detail requires ad- 
dressing a number of differences from the case of normal fluids. In granular fluids, there 
is no "approach to equilibrium" in the usual sense because there is no equilibrium Gibbs 
state, due to the continual loss of energy by inelastic collisions. Thus Eq. (1) constitutes 
a reformulation of Onsager's observation that linear non-equilibrium regression laws can be 
studied via equilibrium fluctuations [14, 18]. It appears that there is a "universal" homo- 
geneous state for an isolated granular fluid replacing the Gibbs state, which is "normal" in 
the sense that all time dependence occurs through the average energy (or, equivalently, the 
granular temperature). This is usually referred to as the Homogeneous Cooling State (HCS), 
and it is the reference state about which spatial perturbations are studied. However, the 
HCS distribution is not simply a function of the global invariants nor is it stationary. The 
time dependence of this reference state, although occurring only through the temperature, 
provides a technical complication in the application of standard methods of non-equilibrium 
statistical mechanics. 

Another related complication is the existence of homogeneous state dynamics (homoge- 
neous perturbations of the HCS) for a granular fluid, even at the hydrodynamic level. In a 
normal fluid, the homogeneous hydrodynamic state is time independent, corresponding to 
the global conserved number, energy, and momentum. These dynamical invariants are the 
hydrodynamic modes in the extreme long wavelength limit. The only dynamical response 
to initial homogeneous perturbations of the equilibrium state corresponds to "microscopic" 
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excitations that decay quickly. Such transient effects can be eliminated from the response 
function in Eq. (1) by a suitable choice for the initial perturbation, such that only the in- 
variants occur in the long wavelength limit. For normal fluids, this is accomplished in all 
previous derivations by choosing an initial local equilibrium ensemble. This choice is both 
practical (eliminating initial homogeneous transients) and physically meaningful. The latter 
refers to the expectation that each cell in a fluid rapidly approaches the equilibrium state for 
that cell, but at its own local density n(r,t), temperature T(r,t), and flow velocity U(r,t). 
Thus most initial preparations will quickly approach the local equilibrium state, and its 
choice as an initial condition avoids the details of these short time transients. A similar 
expectation seems reasonable for granular fluids, where there should be a rapid approach to 
a local HCS. Initial preparation in this state again avoids the problem of initial transients, 
but there is still a residual homogeneous hydrodynamics associated with a uniform pertur- 
bation of all cells. Such excitations correspond to the invariants of a normal fluid, and again 
represent the hydrodynamic modes in the long wavelength limit. They are discussed in the 
next section at the level of phenomenological hydrodynamics. In Sec. Ill, these modes are 
identified in an exact solution to the microscopic Liouville equation, showing the existence of 
hydrodynamic modes in the long wavelength limit and providing the basis for an appropriate 
analysis at finite wavelengths. 

With the homogeneous hydrodynamics characterized, the residual dynamics of the re- 
sponse functions determines the parameters representing effects of spatial gradients. For 
linear response, it is sufficient to consider a single Fourier mode with wavevector of magni- 
tude k = 2tt/\, where A is the wavelength of the spatial perturbation. A formal generator 
for the dynamics of the response function is defined and expanded for long wavelengths to 
second order in k. If the coefficients in this expansion have a finite limit for long times, they 
define the corresponding transport matrix for the phenomenological equations. Accordingly, 
the parameters (cooling rate, pressure, and seven transport coefficients) are given defini- 
tions in terms of correlation functions for the reference HCS in this limit. The analysis is 
straightforward but complex, so many of the details are deferred to appendices. It is perhaps 
helpful to see an overview of the results here before getting too immersed in that analysis 
and associated notation. 

First, the reference HCS is discussed and defined in a representation where it is a station- 
ary solution to a modified Liouville equation, analogous to the equilibrium Gibbs stationary 
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solution. The explicit construction of this solution is a difficult many-body problem and 
is not discussed here beyond noting the many studies of this state by molecular dynamics 
simulation [19]. An exact solution to this Liouville equation for a special homogeneous per- 
turbation of the HCS is constructed and shown to have the same dynamics as that from 
hydrodynamics in the long wavelength limit. This allows identification of the microscopic 
hydrodynamic modes and "invariants" for a granular fluid. A formal solution to the Liou- 
ville equation for corresponding perturbations at finite wavevector is constructed, and the 
response functions defined. The expansion in k of the above solution is performed, with co- 
efficients identified in terms of time correlation functions for the HCS. The phase functions 
defining them are fluxes of the usual conserved densities, and conjugate fluxes associated 
with densities of the invariants. In addition, there are correlation functions for the energy 
loss function due to the inelastic collisions. 

The cooling rate and the pressure in the linear hydrodynamic equations are identified as 
specific averages over the HCS solution. In particular, the pressure is the same average of the 
trace of the microscopic stress tensor as for an equilibrium fluid, but with the equilibrium 
Gibbs distribution replaced by the HCS. This determines the dependence of the pressure on 
density and temperature. The transport coefficients are of two types, those associated with 
the heat and momentum fluxes, and those associated with the cooling rate. In each case 
they can be displayed collectively in matrix form. For the fluxes, they have a Green-Kubo 
representation in terms of the above-mentioned fluxes, 



The first term on the right hand side is a time independent correlation function for the 
HCS. It vanishes for normal fluids with continuous potentials of interaction, and occurs here 
due to the inelasticity of the collisions. Such a term can occur even in the elastic limit for 
singular forces, such as hard spheres [20]. The limit indicated in the second term is the 
usual thermodynamic limit of large volume V and particle number N, followed by the limit 
of large time. The integrand G a p{n,T;i) is a flux-flux correlation function 



The first flux $ a , is one of those associated with the usual densities of number, energy, 
and momentum. The second flux T\, is one of those associated with the densities for the 
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new invariants. Both kind of fluxes are functionals of the phase point T. The evolution 
operator for the dynamics U a p (t, T) is the usual Liouville dynamics, but now compensated 
for the homogeneous dynamics of collisional cooling and homogeneous perturbation. It has 
an invariant subspace {^/ a (T;n,T)}, 

U aP (t, T) ^ (r; n, T) = ^ a (T; n, T) . (4) 

P 

Such a contribution would invalidate the limit in Eq. (2), but it is demonstrated that the 
fluxes Ta are orthogonal to this subspace. This corresponds to the "subtracted fluxes" in 
the Green-Kubo expressions for normal fluids [17], which are recovered from these results in 
the elastic limit. The transport coefficients associated with the cooling rate have a similar 
form, with $ Q replaced by a phase function representing the energy loss from the inelastic 
collisions. 

In both cases, the time correlation functions can be expressed exactly as a time derivative. 
Performing the time integral in Eq. (2) leads to an alternative representation. In fact, three 
representations for the transport coefficients are obtained here. They correspond to the 
familiar results for the diffusion coefficient of an impurity in a normal fluid, 

Id 1 If 1 

D = - lim - — (q • <?o(0) c = - lim 2 (*) • Mt)) c = lim 2 J dt> ( v ° ' u °(*'))c > ( 5 ) 

where qo and Vq are the position and velocity of the impurity, d is the dimension of the 
system, and the angular brackets denote equilibrium averages. The first is the Einstein 
relationship, indicating that the mean square displacement grows linearly with time for long 
times. It was subsequently generalized to other transport coefficients by Helfand [16], and it 
is usually referred to as the Einstein-Helfand representation. The other two representations 
are the intermediate Helfand, and Green-Kubo forms, respectively. 

The derivation of these formal results is accomplished with few restrictions on the dy- 
namics in phase space: deterministic, Markovian, and invertible. This allows a wide range of 
inelastic collision rules currently used for granular fluids, from inelastic hard spheres to soft 
viscoelastic potentials. However, to expose further details of the expressions for the trans- 
port coefficients obtained here, a companion paper following this presentation is specialized 
to the simplest case of smooth, inelastic hard spheres with constant coefficient of restitution 
[21, 22]. Further discussion of the formal results is deferred to the final section. 
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II. PHENOMENOLOGICAL HYDRODYNAMICS 



The purpose of this section is two-fold. First, the phenomenological hydro dynamic equa- 
tions are recalled and the unknown parameters (pressure, cooling rate, transport coefficients) 
indicated. A special solution for spatially homogeneous states is obtained, and the hydro- 
dynamic equations are linearized about that state for small spatial perturbations. A micro- 
scopic representation of these equations is the objective of subsequent sections. The second 
purpose is to characterize the dynamics to be expected from solutions to these equations. 
Specifically, two features new to granular fluids are an inherent time dependence of the coef- 
ficients due to the cooling of the reference state, and a non-trivial dynamics associated with 
homogeneous perturbations of the homogeneous state. In the subsequent sections, identi- 
fying the component of linear response associated with only spatial perturbations requires 
taking explicit account of these two types of homogeneous dynamics. 

Hydrodynamics is a closed description for the dynamics of a fluid in terms of the number 
density n (r, t), the energy density e (r, t), and the momentum density g (r, t). The starting 
point for identifying such a description for a granular fluid is the exact macroscopic balance 
equations for these densities. However, as for normal fluids, it is usual to replace the energy 
density and momentum density by the temperature T (r, t) and flow velocity U (r, t) as the 
hydrodynamic variables, together with the number density. This is accomplished through 
the definitions 



Here m is the mass of a particle and e (n, T) is some specified function of n and T. The 
two most common choices are e (n, T) = dnT/2, where d is the dimension of the system, or 
e (n,T) = e e (n, T), the thermodynamic function for the corresponding equilibrium fluid. 
The former is common in applications of computer simulations (note that the Boltzmann 
constant has been set equal to unity), while the latter is the historical choice in most for- 
mulations of hydrodynamics. Both definitions coincide for the special case of hard spheres. 
For normal and also for granular fluids, the choice made constitutes a definition of tempera- 
ture for non-equilibrium states and has no a priori thermodynamic implications. The exact 
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macroscopic balance equations in terms of n, T, and U are 



An + nV U = 0, (8) 

D t Ui + (mn)- 1 T^Pij = 0, (9) 

3 3 

(iO.(A + o r+ 

where D t = d/dt+XJ- V is the material derivative, C is the cooling rate due to the energy loss 
from the interaction between granular particles, q is the heat flux, and P^ is the pressure 
tensor. These equations have the same form as those for a normal fluid, except for the 
presence of the term involving the cooling rate ( in the equation for the temperature. 

The above balance equations are not a closed set of equations until q, P^, and ( are 
specified as functionals of the hydrodynamic fields, i.e., their space and time dependence 
occurs entirely through these fields. This happens for normal fluids on length and time 
scales long compared to the mean free path and mean free time, respectively, and similar 
conditions may be assumed for granular fluids as well. If, furthermore, the state of the 
system is such that the spatial variation of the fields is smooth, then an expansion of these 
functionals in gradients of the fields can be performed. From fluid symmetry, the results to 
first order in the gradients must have the form: 

Pij = p(n, T)5 l3 - 77(71, T) ^ + ^ - Ik V • Uj - «(n, T)5 tJ V ■ U + . . . , (11) 

q = -X(n,T)VT- fi(n,T)Vn + ..., (12) 

C = Co in, T) + ( u (n, T)V-U + ( T (n, T) V 2 T + C>, T)V 2 n + . . . , (13) 

where 5ij is the Kronecker delta symbol and the dots denote terms proportional to higher 
powers or higher degree spatial derivatives of the hydrodynamic fields than those written 
explicitly. These expressions represent the "constitutive equations" which, together with 
the macroscopic balance equations, give Navier-Stokes order hydrodynamics for a granular 
fluid. Note that the cooling rate is required to second order in the gradients, while the 
pressure tensor and heat flux are required only to first order. The pressure tensor has the 
same form as Newton's viscosity law for a normal fluid, while the expression for the heat 
flux is a generalization of Fourier's law [23, 24]. 



e -n 



de 
dn 



T. 



dUi 



v-c/ + EE^ + v. q = o, 



(10) 
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The expressions (11)-(13) include unspecified functions: the pressure p(n,T), the zeroth 
order in the gradients cooling rate Co (n,T), the transport coefficients from the cooling rate 
( u (n,T), ( T (n,T), ( n (n,T), the shear viscosity rj(n,T), the bulk viscosity k{u,T), the 
thermal conductivity X(n, T), and a new heat flux coefficient /i(n, T). All of these must be 
provided by experiment or by the theoretical justification of this phenomenology. 

Although the Navier-Stokes equations are based on the small gradient forms for the 
constitutive equations, it does not mean that they are limited to systems close to a global 
homogeneous state. Thus, they can be applicable locally over domains larger than the mean 
free path, but the hydrodynamic fields may still vary significantly throughout the system. 
Consequently, a wide range of experimental and simulation conditions for granular fluids 
have been well-described by the Navier-Stokes equations [25-28]. 

The spatially homogeneous solution to Eqs. (8)-(10) for an isolated system (e.g., periodic 
boundary conditions) is 



For simplicity, it will be often considered that U h = in the following, something that 
can always been achieved by means of the appropriate Galilean transformation. All time 
dependence of this state occurs through the homogeneous temperature Th (t) , so this is 
referred to as the homogeneous cooling state (HCS) [29]. Once the functional form for 
Co [nh, T h (t)} has been specified, the first order nonlinear equation (15) can be solved by 
direct integration for a given initial condition. 

Now consider small spatial perturbations of the HCS, assuming that T h (t) is known, 



with 5yp (r, t) sufficiently small that nonlinear terms in the hydrodynamic equations can be 
neglected. The resulting linear equations have coefficients independent of r so the differential 
equations can be given an algebraic representation by means of the Fourier transformation, 



n(r,t) = n h , U(r,t) = U h , T (r, t) — T h (t), 



(14) 



where T h (t) is the solution to 




(15) 



VP ( r > *) = Vp,h + hp (r, t) , {yp} = {n, T, U} , 



(16) 




(17) 
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Moreover, the components of the flow velocity are separated into a longitudinal component 
relative to k, 5U\\ = k ■ 5U, and d— 1 transverse components SU± ji = • 5U, where k = k/k 
and {ef, i — 1, . . . , d — 1} are a set of <i pairwise orthogonal unit vectors. Therefore, from 
now on it is 

{yp} = {n,f,U h U ± }, (18) 

with 5U± denoting the set of the d — 1 components 5U±j. The linearized hydrodynamic 
equations then have the form 

{d t + K! iyd [n h ,T h (t); k]}5y(k,t) = 0, (19) 

where a d + 2 dimensional matrix notation has been introduced for simplicity and d t = d/dt. 
The transport matrix IC hyd is identified as being block diagonal, with decoupled longitudinal 
and transversal submatrices, 

. . / Kll yd \ 
K hyd = 1 , 20 

y K^ d ) 

( n n -in h k \ 



dn h 



- ™) * 2 ^ + - C T T h ) k* -i [c U T h + ^) k 



ip n k iprk 1 



J 



k 2 

(21) 
n h m 

In the above expressions, / is the unit matrix of dimension d — 1. Moreover, /i = e Q + p and 
the following short notations have been introduced: 

'deo\ _ f deo\ _ / <9p\ _ / d£ 



e0 -"l^J T ' e °' T "l5Tj n ' Pn "l^J T ' PT =W) n (23) 
It is understood in Eqs. (21) and (22) that all the quantities are to be evaluated at n = n h 

and T = T h (t). 

There are d+2 independent solutions of Eqs. (19) in terms of which the general solution 
to the initial value problem can be expressed. These are the hydrodynamic modes. For 
normal fluids, they are simply the eigenfunctions of the matrix K, hyd and the hydrodynamic 
excitations are exponential in time with the corresponding eigenvalues. Here, the identifica- 
tion is somewhat less direct due to the dependence of K, hyd on T h (t). For example, the shear 
modes are proportional to 

expi- P4) 
I J n h m J 
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and their time dependence is no longer exponential. 

The general solution to the initial value problem can be written in terms of a matrix of 
response functions (Green's function matrix) corresponding to the representation in Eq. (1), 

5y (fe, t) = C hyd [n h , T h (t); k, t] 5y (fc, 0) , (25) 

which is determined here from the linearized hydrodynamic equations, 

{d t + K. hyd [n h ,T h (t); k]} C hyd [n h ,T h (t); k,t] = 0, (26) 

dff[n h ,T h (0);k,0] = 5 afi . (27) 

Conversely, if C hyd [n h ,T h (t); k,t] were given, the parameters of the transport matrix 
jQhyd [ n/i; T h (t); k] could be determined from 

K hvd [n h ,T h (t); k] = - [d t C hyd [n h ,T h (t); k,t]} C hyd ~ x [n h ,T h (t); k,t] . (28) 

This exposes the formally exact approach to be developed in the next sections to identify 
the parameters of K hyd [n h ,T h (t); k]. First, an exact response function C is defined in place 
of C hyd for 5y (k, t). Next, a transport matrix is defined as in Eq. (28) with C hyd replaced by 
C. On the long time and small k scales for hydrodynamics, this expression must agree with 
JC hyd [rih,Th(i);k], providing a formal definition of the hydrodynamic parameters. Indeed, 
for normal fluids this procedure gives the familiar Helfand and Green-Kubo expressions for 
the transport coefficients in terms of time correlation functions in the equilibrium reference 
state. 

There are two immediate complications for direct extension of this approach to granular 
fluids. The first is the parametrization of IC hyd [n h) T h (t) \ k] by the time dependent temper- 
ature Th(t). For example, the transport coefficients are functions of this temperature. Since 
Th(t) is determined autonomously from Eq. (15), it is useful to suppress this dynamics by 
the definitions 



JC hyd [n h ,T h (t);k] = [)C hyd (n,T;k)} 



n=n h ,T=T h (t) 



C h y d [n h ,T h (t);k,t] = C hyd (n,T;k,t) . (29) 

L 1 n=n h ,T=T h (t) 

Thus IC hyd (n, T; k) and C hyd (n, T; k, t) are the universal function associated with a general 
class of reference states. However, in this form t and T become independent variables so, 
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for example, the equation for C hyd (n, T; k, t) is not Eq. (26) but 



d t -( (n,T)T-^ + IC hyd (n,T;k) 



C hyd (n, T; k, t) = 0, 



(30) 



again with the initial condition C^jf (n, T; k, 0) = 5 a p. The time derivative is now at constant 
T. The new term with the temperature derivative in Eq. (30) is the generator for the 
dynamics of T h (t) in the sense that for any arbitrary function X(T) it is 



X[T h {t-T)\ =exp 



-t( (n h ,T)T^ 



X(T), 



(31) 



where T h (t;T) is the solution to Eq. (15) with 7^(0) = T. The proof is given in Appendix 
A. Equation (30) now has time independent coefficients, but at the price of introducing 
a new independent variable T. Still, it is most convenient to work with the generic forms 
IC hyd (n, T; k) and C hyd (n, T; k, t) and the counterpart of Eq. (28), 



IC hyd (n, T; k) = 



dt - Co (n, T) T-^ 



C hyd (n, T; k, t) \ C hyd ~ l (n, T; k, t) . (32) 



The second complication is the existence of a non-zero generator for homogeneous (zero 
wave vector) dynamics, i.e. K hyd (n,T;0) ^ 0, or C^f (n,T;0,t) ^ 5 afS . The latter can be 
calculated directly to give (see Appendix A): 



C hyd (n,T;0,t) 



C hyd Q 
/ 



(33) 



where 



C} yd (n,T; 0,t) 



1 



\dn)T(-t 





dT 



9T \ 








\ 



V 



(34) 



/ 



] T(-t;T) \dT(- 

I 

and I is again the unit matrix of dimension d — 1. The interpretation of the above result is 
clear when evaluated at n = rih, T = Th(t), where the non trivial matrix elements become 



C% d [n h ,T h (t);0,t] = 



dn h 



T h (0) 



C% d [n h ,T h (t);0,t] = 



(35) 



This is just the linear response of the solution to Eq. (15) due to changes in the initial 
conditions. It is the additional dynamics of the temperature beyond that of T h (t), due 
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to initial homogeneous density and temperature perturbations, as the system attempts to 
approach a new HCS. Note that this it is a hydrodynamic excitation, and is distinct from 
the rapid homogeneous relaxation of other microscopic modes on a much shorter time scale 
outside the scope of hydrodynamics. Hence for granular fluids the hydrodynamic modes 
cannot be identified simply as those that vanish for k — > 0, as for a normal fluid. Instead, 
they are identified as those modes whose homogeneous dynamics becomes that of Eq. (33) 
for k — > 0. This is done in the next section. 

III. STATISTICAL MECHANICS 

Consider a volume V enclosing N particles that interact and they lose energy as a result 
of this interaction. Also suppose that the interactions are specified in such a way that if 
the positions and velocities of each of the particles are given at a time to, then there exists 
a well defined trajectory for the evolution of the system for all times both earlier and later 
than t - The microscopic initial state of the system is entirely determined by the positions 
and velocities of all particles, denoted by a point in a 2dN dimensional phase space T = 
{q r , v r ; r — 1, . . . , N}. The state of the system at a later time t is completely characterized 
by the positions and velocities of all particles at that time r t = {q r (t),v r (t);r = 1, . . . , iV}. 
In summary, the dynamics is Markovian and invertible. 

The statistical mechanics for this system is comprised of the dynamics just described, a 
macrostate specified in terms of a probability density p(T), and a set of observables (mea- 
surables). The expectation value for an observable A at time t > 0, given a state p(T) at 
t — 0, is defined by 



The notation A(T t ) indicates the function of a given phase point shifted forward in time along 
a trajectory. Equivalently, this can be considered as a function of the initial phase point 
and the time, A(T,t). This second representation allows the introduction of a generator L 
for the time dependence defined by 



where the explicit form of L is determined by the specific rule chosen for the interactions 
among the particles. This is left unspecified at the moment. Due to the assumption of 




(36) 




(37) 
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invertibility made above, the dynamics can be transferred from the observable A(T) to the 
state p(r), by defining an adjoint generator L, 

J dT p{T)e tL A{T) = J dT [e-* r p(r)j A(F) = J dT p(T, t)A(T). (38) 

The form generated by L is referred to as Liouville dynamics. This equivalence of observable 
and state dynamics is expressed in the above notation as 

(A(t);0) = (A;t}. (39) 

In summary, the dynamics of phase functions is governed by an equation of the form 

(d t -L)A(T,t) = 0, (40) 
and that of the probability distribution in phase space by a Liouville equation 

(d t + L)p(T,t)=0, (41) 

where L and L are time independent operators. 

Time correlation functions are defined in a similar way: 

Cab (t) = (A(t)B(p); 0) ee J dT p(T)A(T t )B(T). (42) 

In terms of the generators introduced above, the correlation functions can be written 

C AB (t) = J dTp(T) [e tL A(T)] B(T), (43) 

or, equivalently, 

C AB (t) = J dTA(T)e- tT [p(T)B(T)] . (44) 
Further comment on these generators and some examples are given in Appendix B. 

A. Homogeneous Reference State 

In contrast to normal fluids, there is no stationary solution to the Liouville equation (41) 
for an isolated granular fluid, because the average energy of the system decreases with time, 

d t (E-t) = (LE;t}<0, (45) 
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where E(T) is the phase function corresponding to the total energy and the equal sign cor- 
responds to the elastic limit. This loss of energy due to interactions among the particles is a 
primary characteristic of granular fluids. It is convenient to introduce a granular tempera- 
ture instead of the average energy in the same way as is done in (6) for the phenomenological 
equations. For a homogeneous state, the flow velocity U is uniform and can be chosen to 
vanish by an appropriate Galilean transformation, as already mentioned. The temperature 
definition for the homogeneous state then is 

e [n,T(t)] = ^{E;t), (46) 

and Eq. (45) becomes 

d t T(t) = -((t)T(t), (47) 
where the "cooling rate" is identified as 

f^-MSj iS;i>> -°' (48) 

Upon writing the last inequality, it has been taken into account that for both choices 
of the function e (n, T) discussed in the previous section, e is an increasing function of T. 
The above shows that there is no "approach to equilibrium" for a granular fluid, except in 
the elastic collision limit where ( (t) = 0, as no such stationary state exists for an isolated 
system. Consequently, there is a large class of dynamical, homogeneous states depending 
on the initial preparation. However, as in a normal fluid, it is expected (and observed in 
computer simulations [19]) that there is a rapid relaxation of velocities after a few collisions 
to a "universal" state whose entire time dependence occurs through the cooling temperature. 
This special state is called the homogeneous cooling state (HCS), since it is the microscopic 
analog of the hydrodynamic HCS described in Sec. II above. The probability density for 
the HCS then becomes [3, 4] 

p h (T,t) = p h [T;n h ,T h (t)}, (49) 

where rih is the uniform, time independent density and T h (t) the decreasing temperature, 
consistently with the notation used in Sec. II. Because all time dependence occurs through 
T h (t), it follows that ( h (t) = ( [n h ,T h (t)}, whose form is obtained from 

Co (n, T) = --L J dV Ph (T; n, T)LE (T) . (50) 
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The functional form of ph(T; n, T, U) is determined self-consistently by the Liouville equation 
(41) or, equivalently, using the results in Appendix A, 

C TPh (T;n,T,U)=0, (51) 

with the definition 

Z T = -Co (n, T) T A + L. (52) 

Here, a homogeneous velocity field [/ is formally considered, since it is convenient for 
later purposes. Its presence does not affect the value of the integral on the right hand 
side of Eq. (50). Note that since the HCS is normal, i.e. all the time dependence occurs 
through the temperature, the time derivative is replaced by the scaling operator for cooling, 
—Co (n,T)Td/dT, of Eq. (31). Equations (50) and (51) constitute a pair of time indepen- 
dent equations to determine both ph(X; n, T, U) and ( (n,T). Once the latter is known, 
T h (t) is determined from the solution to Eq. (47) which becomes 

d t T h (t) = -( [n h ,T h (t)]T h (t). (53) 

Finally, Ph{F; n, T, U) is also evaluated at n = nh, T = T^t), and U = Uh = 0. 

This definition of the HCS in terms of the solution to the Liouville equation provides the 
first of the unknown hydrodynamic parameters, the zeroth order cooling rate Co [nh,Th(t)]. 
It now has a precise and exact definition, Eq. (50), from which to determine its density and 
temperature dependence. 

It is easily verified that the equilibrium probability densities for normal fluids (e.g., Gibbs 
ensembles or maximum entropy ensembles for conserved densities) are not solutions to Eq. 
(51). Thus, the effects of inelastic collisions are two-fold. First, they introduce an inherent 
time dependence due to energy loss, through Th(t), and second, they change the form of 
the probability density in a way that cannot be simply represented by the global invariants 
for a normal fluid. For the analysis here, it is assumed that Eq. (51) can be solved to good 
approximation and its solution is taken as known. The existence of this solution is supported 
by molecular dynamics simulations that show a rapid approach to a state whose granular 
temperature obeys Eq. (53), with uniform density and flow velocity [19]. Furthermore, at 
sufficiently low density, the reduced distribution function obtained by integrating out all de- 
grees of freedom but one in Eqs. (50) and (51), should be approximated by the corresponding 
Boltzmann limit. This limit also supports a HCS solution, as verified by simulations of the 
Boltzmann equation using the direct simulation Monte Carlo method [30]. 
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B. A Class of Homogeneous Solutions to the Liouville Equation 



Since cooling is inherent to the inelastic collisions for all states, it is useful to consider 
the time dependence of a general state comprised of that through T h (t) plus any residual 
time dependence, i.e., p (T; t) = p[T;n h ,T h (t),U;t\. Then the Liouville equation (41) takes 
the form 

(d t + C T )p(T;n,T,U;t) = 0. (54) 

It is again understood that t and T are now independent variables and the time derivative is 
taken at constant T. The new generator for the dynamics Ct, given in Eq. (52), incorporates 
both that for the trajectories in phase space and the effects of cooling. The final solution 
to this equation is to be evaluated at n — rih, T = Th(t). In this representation, the 
distribution function of the HCS, ph(T;n,T,U), is seen to be a stationary solution to the 
Liouville equation, as shown by Eq. (51). 

Once p h (T;n,T, U) is known, a class of other solutions can be constructed for homoge- 
neous perturbations of the HCS. These are the microscopic analogues of the hydrodynamic 
solution defined by C hyd (n, T; 0, t) given in Eq. (33). To see this consider the initial condition 

p(T;0) = Ph (T-n + 5n,T + 6T,U + SU) 

~ Ph (r; n, T,U) + J2 (F; n, T, U) 6y a (0) , (55) 

a 

with the spatially uniform perturbations {5y a (0)} = {5n,5T,5U} and 

iTf (r T TT\ ( d P h ( r ' n > T > U ) \ fti R\ 

^(r;n,T,C/)=^ — j . (56) 

The corresponding solution to the Liouville equation (54) is 



p(T,t) = e- CTt 



Ph 



(r ; n, t,u) + J2 ( r ; ^ T > u) s Va (o) 



= Ph (r ; n, T, U) + e~ CTt (r; n, T, U) 5y a (0) . (57) 

a 

The second equality follows from the stationarity of ph for the generator Ct, Eq. (51). The 
second term in the last expression can be evaluated using the identity 

C T ^ a (r; n, T,U) = J2 */3 ( r ; ^ T > U) /Cg d (n, T; 0), (58) 
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which follows by direct calculation using the definition of Ct and Here IC hyd (n,T; 0) is 
the same matrix as obtained by evaluating Eqs. (20)-(22) at k = 0. The proof is given in 
Appendix C, where it is also shown that Eq. (58) in turn gives 

e~ Z ^ a (r; n, T, U) = ^ 9p (r ; n, T, C7) (n, T; 0, f ) . (59) 

P 

Here (n, T; 0, t) is the same hydrodynamic response matrix of Eq. (33). The exact 
solution to the Liouville equation is therefore 

p (r; f) = Ph (r; n, T, C7) + ^ ^ a (r; n, T, £/) fy a (t) , (60) 

a 

where 

^(t) = ^^(n,T;0,t)^(0). (61) 

The special choice of the \& a for initial perturbations is seen to excite only hydrodynamic 
modes, and no other microscopic homogeneous excitations. In this sense, the ^ a (I"; n, T, U) 
can be considered the microscopic hydrodynamic modes in the long wavelength limit. For a 
normal fluid, they become functions of the global invariants and therefore time independent, 
which are indeed the hydrodynamic excitation at k = in that case. For a granular fluid, as 
discussed in the previous section, there is a nontrivial dynamics even at k — 0. The analogy 
can be made more direct by rewriting Eq. (59) in the form 

E u ^ (*> T ) *f> ( r ; n > T > u ) = ( r 5 n ' T ' u ) ' ( 62 ) 

where the new matrix evolution operator U (t, T) is defined as 

U aP (t, T) = Cjf- 1 (n, T; 0, t) e"^'. (63) 

The dynamics described by U (t, T) is the Liouville dynamics, compensated for both effects 
of cooling, through the scaling generator in £ T , and the dynamics of homogeneous pertur- 
bations, through the response function C hyd (n, T; 0,t). Consequently, U (t,T) provides the 
dynamics associated with spatial perturbations. It will be seen below that this operator de- 
fines the time dependence of the correlation functions representing all transport coefficients. 
Equation (62) shows that the functions {^ a } are the global invariants for this dynamics. 
Note that although it is convenient to keep U ^ at a formal level for the discussion here, 
there is no problem in taking the limit U — > in all the obtained results. 
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IV. INITIAL PREPARATION AND LINEAR RESPONSE 



In this section, the response of the system to an initial spatial perturbation in the hydro- 
dynamic fields relative to the HCS is studied, in order to extract the exact analogue of the 
hydrodynamic transport matrix described on phenomenological grounds in Eqs. (20)-(22) 
above. Consider an initial perturbation generalizing (55) to inhomogeneous states, 



p (r; 0) = Ph (r; n,T) + J2jdr ^ (r ; n, T; r) 8y a (r, 0) , 



(64) 



where now {5y a (r,0)} = {8n(r, 0), 8T{r, 0), 8U(r, 0)} are the initial space dependent de- 
viations of the hydrodynamic fields from their values in the reference HCS and we have 
taken for simplicity U — 0. Thus the argument U will be omitted in the following when it 
vanishes. The {ip a } are spatial densities corresponding to the invariants {^ a }, 



V a (T; n, T) = J dr^ a (F;n,T;r). 



(65) 



This choice of this initial condition assures that the long wavelength limit of the perturbation 
gives the hydrodynamic solution (60). The hydrodynamic fields {8y a (r,t)} as identified 
below Eq. (64), are averages of corresponding phase functions {a a (T; n, T; r)}, or, more 
precisely, it is 



Sy a (r, t) = J dT [p (r; t) - p h (T; n, T)] a a (T; n, T; r) , 
K (r; n, T;r)}= [jsf (r ; r), — [£ (r ; r) - e ^ (r ; r)] 6[l: ' ] 



nm 



(66) 



(67) 



The local microscopic phase functions for the number density, N "(r), energy density, £ (r), 
and momentum density, C/ (r), are 

/ i \ 



^ AT(r;r)^ 



iV 



^ 5 (r - q r ) 



r=l 



"2" 



+ | Er^s V (Qrs) 



(68) 



/ 



In this expression, q rs = q r — q s and V (q rs ) is the pair potential for the conservative part of 
the interaction among particles, as discussed in Appendix B. Normalization of both p (T; 0) 
and ph (T; n, T), and the representation (66) for 8y a (r, 0) require 

J dT i, a (r; n, T; r) = 0, ^ rfr a a (r ; n, T; r) ^ (r ; n, T; r') = 5^5 (r - r') . (69) 
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The second relation above shows that {ap} and {ipp} comprise two biorthogonal sets. 

The conditions (65) and (69) are satisfied if the densities ip a (T; r) are generated from a 
local HCS distribution, p ih (T \ {yp}), through 

'Sm (r\{yp})~ 



^ a (T;n,T; r) 



(70) 

{yp}={n,T,0} 



$Va (r) 

This local HCS distribution is the analogue of the local equilibrium distribution for a normal 
fluid. Qualitatively, it corresponds to the condition that each local cell has an HCS distri- 
bution characterized by the parameters {y a (r, 0) = y a ^ + fya ( r , 0)}. This characterization 
as a local form for the HCS requires 



(r, 0) = J alT [p lh (r | {^(0)}) - p h (F, n, T)] a a (T; {y p } ; r) , (71) 



h-I 



Pih (r| {y a ,h}) = Ph (r; {y a ,h}) , (72) 



dr r , 



S n Pih 



d n p h (T;n,T, U) 



(73) 



5y a (n) ...6yp (r n ) J {y0}={n<TjU} dy a>h . . . dy^ h 
The first equality states that the local state has the exact average values for the {a a }; the 
other two refer to the uniform limit. This is sufficient for the conditions (65) and (69) to be 
verified. 

The formal solution to the Liouville equation for this initial condition, generalizing Eq. 
(57) to spatially varying perturbations, is 

p (r; t) = p h (r; n,T) + J2 fdr e~ Z ^ a (T; n, T; r) 5y a (r, 0) . (74) 
The response in the hydrodynamic fields is therefore, in matrix notation, 

5y (r, t)= J dr' C (n, T;r- r', t) 6y (r', 0) , (75) 
with the matrix of response functions defined by 

C af3 (n, T;r- r\ t) = j dT a a (T; n, T; r) e"^*^ (r; n, T, r') . (76) 

This notation expresses the translational invariance of the reference HCS and the generator 
for the dynamics. This is the exact response function for the chosen initial perturbation, 
representing all possible excitations of the many-body dynamics. The corresponding Fourier 
representation is 

C aP (n, T; fe, t) = V- 1 J dTa a (r ; n, T; k) e'^MT; n, T, -k). (77) 
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These response functions will be the primary objects of study in all of the following. 
As a consequence of Eq. (65), at k = the ^(r; n, T, fe)'s become the invariants 



: (r; n, T; 0) = J dr ^ a (T; n, T; r) = V a (T; n, T) . 



(78) 



Therefore, it follows from Eqs. (59) and (69) that the exact response function in the long 
wavelength limit is the same as that from hydrodynamics, 



C aP (n,T; Q,t) = C*f{n,T-Q,t), 



(79) 



with C^ d (n,T;0,t) given by Eq. (33). By construction, the choice of initial preparation 
made has eliminated all microscopic homogeneous transients and only the hydrodynamic 
mode is excited. This result is exact, and shows that the microscopic dynamics supports a 
hydrodynamic response at long wavelengths. At this formal level, analyticity in k is sufficient 
to admit the possibility of Navier-Stokes modes. 

To identify the macroscopic hydrodynamic equations, it is useful first to write an exact 
equation for the response function in a form similar to (30), 



d t -( (n,T)T-^+lC(n,T;k,t) 



C (n, T; k, t) = 0, C a p (n, T; k, 0) = 5 a p, (80) 



which provides a formal expression for the generalized transport matrix, 

d 



K{n,T\k,t) = - 



d t -( (n,T)T 



K hyd {n, T; 0) 



dT 



d 



C(n,T;k,t) > C _1 (n, T; k, t) 



d t - Co (n, T) T— + K** (n, T; 0) 



C(n,T;k,t) } C' 1 (n,T;k,t). 

(81) 



The contribution from k = has been extracted explicitly, since Eq. (79) implies it is exact at 
all times. The other term is proportional to the time derivative relevant for the homogeneous 
dynamics and therefore is of order k. The full hydrodynamic matrix, as given by Eqs. (20)- 
(22), when it exists, follows from this formal result for small k (long wavelengths) and long 
times, 

K, hyd (n, T; k) = lim /C (n, T; k, t) . (82) 

t»t ,k«k 

The characteristic time t and wavelength A;^ 1 are expected to be the mean free time and 
mean free path, respectively. Comparison of this expression with the form (20) provides a 
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"derivation" of the linear hydrodynamic equations, and also gives the coefficients of those 
equations in terms of the response functions. A detailed comparison up through order k 2 is 
the objective of the next few sections. 

The result (81) is the first of three exact representations of the transport matrix to be 
obtained here. Its expansion to order k 2 leads directly to the Einstein-Helfand representa- 
tion of the transport coefficients as the long time limit of time derivatives of correlation 
functions. This is analogous to the diffusion coefficient D represented in terms of the 
time derivative of the mean square displacement, i.e. the first representation in Eq. (5) 
[31]. However, due to the homogeneous state dynamics, the relevant time derivative is 
[d t - Co (n, T) Td/dT + K hyd (n, T; 0)] , so this form may not be optimal in practice. There 
is an intermediate Helfand form which entails correlation functions with non-zero long time 
limits determining the transport coefficients. Finally, the third equivalent representation is 
the Green-Kubo form in terms of time integrals of correlation functions. These second and 
third forms are given in the next section and utilized to implement the k expansion. 

V. NAVIER-STOKES HYDRODYNAMICS 

The linearized Navier-Stokes equations follow from an evaluation of the transport matrix 
of Eq. (81) to order k 2 . This can be accomplished by a direct expansion of C Q/ 3 (n, T; k, t) in 
powers of k [31]. It is somewhat more instructive to proceed in a different manner, using the 
microscopic conservation laws to expose the dominant k dependence. This allows interpre- 
tation of the phase functions occurring in the correlation functions of the final expressions. 

A. Consequences of Conservation Laws 

For normal fluids, the variables a a (T;n,T;k) are the Fourier transforms of linear 
combinations of the local conserved densities, so their time derivatives are equal to 
ik ■ f a (T;n,T;k), where the f a are the associated microscopic fluxes. The proportion- 
ality to k of the time derivatives means that they vanish in the long wavelength limit, as 
appropriate for a conserved density. This allows evaluation of the time derivative in Eq. 
(81) and shows the transport matrix is of order k. Then shifting the time dependence to the 
other density in the response function and using again the conservation law, the dependence 
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through order k 2 is exposed in terms of correlation functions involving the fluxes [17]. It is 
somewhat more complicated for granular fluids, although the general idea is the same. 

Consider first the time derivative occurring in Eq. (81). Using the definition of the matrix 
of response functions, Eq. (77), it can be transformed into 



--7 



dT 



d t -( (n,T)T-^ + JC h y d (n,T;0) 



d 



C(n,T;k,t) 



L - Co (n, T) T— + /C^ (n, T; 0) 



a (T; n, T; k)j e -£T fy(r; n, T; -k), 

(83) 

where the adjoint generator L of L has been introduced. As mentioned above, for normal 
fluids it is La a — ik • f a (T;n,T; k). Here, the inelastic collisions give rise to an additional 
energy loss w(T; k), which is not proportional to k and therefore cannot be absorbed in the 
flux. The new relationships are (see Appendix D): 

w(T;k) 



La a (T; n, T; k) = ik ■ f a (T; n, T; k) - 5 a2 



eo,r 



(84) 



The detailed expressions of f a (T; n, T; k) and w(T; k) are given in Appendix D. The pre- 
factor l/e ,T m the energy loss term, appears because of the definition of a 2 . For a ^ 2, 
the right side gives the usual fluxes for number and momentum density. Inclusion of the 
additional terms -Co (n, T) Td/dT + K, hyd (n, T; 0) in Eq. (83) modifies this result to 

d ~ 



L-( (n,T)T 



dT 



a a (T; n, T; k) + ^ 1C*f (n, T; 0) n, T; k) 



= ik f a (r; n, T; k) - 5 a2 £{T; n, T; k), 
where £(T; n, T; k) is defined by 
1 



£(T;n,T;k) = 



e ,T 



w 



T; k) - V- 1 J2 (r; n, T; k) J dT ^ a (T; n, T) w(T; 0) 



(85) 



(86) 



The first contribution in this expression is the phase function whose average in the HCS 
gives the cooling rate, 



— [dT Ph (r; n, T) w{T- 0) = Co(n, T)T. (87) 
eo,r J 

The remaining terms assure that £(T;n,T; 0) is orthogonal to the invariants, namely that 

(88) 



£(T; n, T; 0) = (l - pt) ^1 = - (l - pt) LE ^ 



eo,T 
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Here, Pt is the projection operator onto the set {a a (T; n, T, 0)}, 

ptX(r) = \/- 1 ^a a (r;n,T;0) f dT^ a (T;n,T) X(T). (89) 

To verify that P^ is really a projection operator, recall that as a consequence of Eq. (69), 
{a a (T; n, T, ; 0)} and {^/ a (T; n, T)} form a biorthogonal set, 

V' 1 J dTa a (T;n,T;0)^^T;n,T) = 5 a/3 . (90) 

Use of Eq. (85) in Eq. (83) shows that the correlation functions C a p{n,T;k,t) obey the 
equations 



d t - CoT A + /C^ d (n, T; 0) 



C (n, T; k,t) = ik-D (n, T; k,t) - S (n, T; k, t) . (91) 



The new correlation functions on the right hand side, D a/3 (n, T; k,t) and S a p(n,T;k,t), 
are similar to C a p (n, T; k, t) but with a a replaced by f a and 5 a 2i, respectively, 

D aP {n, T; k, t) = V- 1 J dT f a (T; n, T; fe^V^r; n, T; -k), (92) 

S aP (n,T;k,t) = 5 a2 V- 1 J dT£(T;n,T;k)e- XTt MT;n,T;-k). (93) 

The utility of Eq. (91) is that its use in Eq. (81) leads to an expression in which the transport 
matrix is exposed exactly through first order in k, 

K {n, T; k, t) = JC hyd (n, T; 0) - \ik ■ D (n, T; k,t)-S (n, T; k, t)j C' 1 (n, T; k, t) . (94) 

It follows from Eqs. (59) and (88) that S (n, T; 0, t) = 0, so the term between square 
brackets on the right hand side of Eq. (94) is at least of order k, as said above. However, this 
representation is still not optimal since the right hand side has the homogeneous dynamics 
of C (n, T; 0, t) that should be cancelled. This technical point is addressed by incorporating 
C^ 1 (n, T; 0, t) in the evolution operator for the correlation functions by using U a p (t, T) 
introduced in Eq. (63), i.e. defining 

MT-n,T;k,t) = J2^(t,T)^(T-n,T;k). (95) 

There is no k = dynamics for ip a (T; n, T; k, t) since ip a (T; n, T; 0) is an invariant; see Eq. 
(78). The transport matrix (94) then becomes 

K(n, T; k,t) = K, hyd {n, T; 0) - [ik ■ D(n, T; k, t) - S(n, T; k, t)] C~\n, T; k, t). (96) 
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The correlation functions with the over-bar are the same as those with the tilde, except that 



now are defined with the dynamics of (95), 

C a p(n, T; k, t) = V- 1 J dTa a (T; n, T; k)MV; n, T; -k, t), (97) 

D a/3 (n, T- k, t) = V- 1 J dT f a (T; n, T; k)^(T; n, T; -k, t), (98) 

S a() (n, T; k, t) = 5 a2 V- 1 [ dT £(T; n, T; k^T; n, T; -k, t). (99) 



Equation (96) gives the intermediate Helfand representation referred to at the end of the 
last section. It has the advantage of being expressed in terms of the appropriate dynamics 
of U a/ 3 (t,T), as well as avoiding the complex time derivative of Eq. (81). 

The equivalent Green-Kubo form is obtained by representing the correlation functions in 
Eq. (96) as time integrals. This is accomplished by observing that there are "conjugate" 
conservation laws associated with ^(r; n, T; k, t). Their existence follows from the fact that 
the ipaS are the densities associated with the invariants ^/ a 's. The conjugate conservation 
laws are 

d t MT; n, T; k, t) - ik ■ % (T; n, T; k, t) = 0. (100) 
The new fluxes 7 a are identified in Appendix D as 

j a (r; n, T; k,t) = J2 U a p (t, T) 7/3 (r ; n, T; k) (101) 



with 

ik ■ j a (r; n, T; k) = -CtMT; n, T; k) + ]T K h £{n, T; 0)MT; n, T; k). (102) 

These new conservation laws give directly 

d t C a/3 (n, T; k, t)+ik- E aP {n, T; k, t) = 0, (103) 

d t D a p(n, T; k, t) + ik ■ F Q/3 (n, T; k, t) = 0, (104) 

d t S a p(n, T; k,t) + ik ■ N aP {n, T; k, t) = 0, (105) 

where 

E af3 (n,T;k,t) = V- 1 [ dra a (F; n, T; fc)7 /3 (r; n, T; —k, t), (106) 
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F Q/3 (n, T; k, t) = V' 1 j dT f a (T; n, T; k)^(T; n, T; -k, t), 
N a/3 (n,T;k,t)=5 a2 V- 1 J dT£(T;n,T;k)j p (T;n,T;-k,t). 



(107) 
(108) 



Note that F a/ g is a second-rank tensor. Integrating Eqs. (103)-(105) allows C(n,T;k,t), 
D(n,T; k,t), and S(n,T;k,t) to be eliminated from Eq. (96) in favor of E, F, and N, 
exposing a higher order dependence on k, 

JC(n,T;k,t) = JC hyd (n,T;0) - [ik -D(n,T;k,0) 

-S{n,T;k,0) + ik- [ dt'N(n,T;k,t') + kk: [ dt'F(n,T;k,t') 
Jo Jo 

x I + ik- dt , (T 1 (n,T\k,lf)E(n,T\k,lf)(r 1 (n,T\kJ) 
Jo 



(109) 



This is the Green-Kubo form for the transport matrix. An advantage of this form is a further 
exposure of the explicit dependence on k. In both Eqs. (96) and (109), relevant correlation 
functions are seen to be those composed from the conserved densities |a Q ,^ a |, the fluxes 
of the two kinds of conservation laws < f a ,la \ , and the source term for inelastic collisions 



£. All of the time dependence is given by the evolution operator U a p (t, T) which is that for 
the N particle motion in phase space, but compensated for all homogeneous dynamics. 



B. Green-Kubo Form to Order k 2 



Retaining only contributions up through order k 2 in Eq. (109) gives 

K(n, T; k, t) = lC hyd (n, T; 0) - ik \k ■ D(n, T; 0, 0) + Z(n, T) 
+k 2 [A(n,T) + Y(n,T)] , 



(110) 



where the meaning of the different terms will be discussed next. The first order in k terms on 
the right hand side of this equation provide the parameters for Euler order hydrodynamics. 
At this order, the susceptibilities (pressure and pressure derivatives) are defined in terms of 
the time independent correlation function D a p{n, T; 0, 0), while the transport coefficient ( u 
is given by the Green-Kubo expression 



Z a0 (n,T)=5 a2 5 f33 T( u (n,T), 



T( u (n,T) = -k 



(n, T; 0) - lim / dt'N 23 (n, T; 0,0 
o 



(111) 
(112) 
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The above identification has been made by comparison of the expression obtained here with 
the phenomenological transport matrix in Eqs. (20)-(22). Here and below, the notation for 
Taylor series expansion of any function X(k) is 



X(k) = X(0) + ik ■ - kk : X (2) + 



(113) 



At order k 2 , the Navier-Stokes transport coefficients in Eq. (110) are of two types. The 
first type are those obtained from A(n, T), and represent the dissipative contributions to 
the fluxes. They correspond to the shear and bulk viscosities, thermal conductivity, and ji 
coefficient in Eq. (12). In the Green-Kubo form they are determined by 

ft 



A Q/3 (n,T) = kk 



D% {n, T; 0) - lim / dt' G Q/3 (n, T; if) 



(114) 



with 



G(n, T; t) = F(n, T; 0, t) - D(n, T; 0, 0)E(n, T; 0, t). (115) 

The transport coefficients of the second kind are those following from Y(n,T) in Eq. (110) 
and represent the second order gradient contributions to the cooling rate, i.e. the coefficients 
C T and C" in Eq. (13), 



Y a(3 (n,T) = -5 a2 kk : 



S ( ^{n,T;0) -lim / dt' H 2j a(n, T; 0, t') 
o 



(116) 



■ (!)/ 



H 2/3 (n, T; 0, t) = W 2 >(n, T; t) + T( u (n, T) kE w (n, T; 0, t) . (117) 

The superscripts (1) and (2) denote coefficients in the expansion of correlation functions in 
powers of ik, as indicated in Eq. (113). 



C. Intermediate Helfand Form to Order k 2 



The intermediate Helfand form to order k 2 follows from direct expansion of Eq. (96) 
The structure is the same as in Eq. (110) as well as the contribution from D a p(n,T; 0,t) = 
D a/3 (n,T; 0, 0) . On the other hand, the transport coefficients are now given by 

T( u (n,T) = -Kmk-S%(n,T;t), 



A(n, T) = limfcfc : [D (1) (n,T;t) - D(n, T; 0, 0)C^\n, T; t) I , 



Y a p(n, T) = -5 a2 lim kk : sg (n, T; f ) + T( u (n, T) k ■ (n, T;t)]. 



(118) 

(119) 
(120) 
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The equivalence of these results with the Green-Kubo forms given in the previous subsection 
can be seen by noting that the conservation laws of Eqs. (103)-(105) to first order in k give 

E(n,T;0,t) = -d t C (1 \n,T;t), F(n, T; 0, t) = —d t D^ (n, T; t) , (121) 

N(n,T;0,t) = -d t S {1 \n,T;t), N (1) (n, T; t) = -d t S (2) {n, T; t). (122) 

These allows the time integrals in the Green-Kubo expressions to be performed, giving 
directly Eqs. (118)-(120). 

D. Einstein-Helfand Form to Order k 2 



Finally, the Einstein-Helfand form to order k 2 follows from direct expansion of Eq. (81), 
K{n, T; k,t)=K, {n, T; 0) + ik ■ /C (1) (n, T; t) - kk : K (2) (n, T; t) , (123) 



with 

7C (1) (n,T;t) = - 



dt - Co (n, T) T-^ + K h * d in, T; 0) 



C (1) i^T-t^C- 1 in,T;0;t) 

(124) 



and 



K( 2 ) in,T;t) 



d 



d t - Co in, T) T-^ + K!* in, T; 0) 



C (2) (n,T;t) 



+A: (1) in, T; t) C (1) in, T; t) } C' 1 in, T, 0; t) . 



(125) 



This form does not separate explicitly the contributions leading to the transport coefficients 
at both Euler and Navier-Stokes orders for the cooling rate. In the previous two represen- 
tations, this was possible because the microscopic phase function £(T; n, T; k) due to energy 
loss in the inelastic collisions appears explicitly. 



E. Dynamics and Projected Fluxes 

The Green-Kubo expressions involve the long time limit of time integrals over correlation 
functions. This presumes the correlation functions decay sufficiently fast for the integrals 
to exist. This decay time sets the time scale after which the hydrodynamic description can 
apply. If these integrals converge on that time scale then the Helfand formulas also reach 
their limiting plateau values on the same time scale. 
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To explore this time dependence further, consider the correlation function characterizing 
the transport coefficients A aj g(n, T) associated with the heat and momentum fluxes (see Eq. 
(H4)), 

G(ra,T;t) = F(n, T; 0, t) — ~D(n, T; 0, 0)E(n, T; 0, t) 
= V- 1 j dTf(T;n,T;0)[j(T;n,T;0,t) 

-i){T- n, T; O)^ 1 J dTa(T; n, T; 0)7(1"; n, T; 0, t) 

= V- 1 J dT/(r; n, T; 0) (1 - P) 7(r; n, T; 0, t). (126) 

Then both contributions to G combine to form the projected part of the fluxes (1 — P) 7, 
where P is the projection onto the set {^ a (T;n,T)}, 

PX(T) = ^ _1 *(r;n,T) J dTa(T;n,T;0)X{T) 

= V- 1 ^«( r ; ^ T ) j dT a Q (r ; n, T; 0)X(r). (127) 
Thus (1 — P) is a projection orthogonal to the invariants, and (1 — P)U (t,T) P = 0, so 

(1 - P) W (f , T) 7(r ; n, T) = (1 - P) U (t, T) (1 - P) 7(1^ n, T). (128) 
The time correlation function (126) then can be rewritten as 

G a/3 (n,T;t) = V- 1 J dT^ a (T;n,T)[U(t,T)T(T;n,T)] p 

= V- 1 I dr *«( r ; ^ T )Upx (t, T) T A (r ; n, T), (129) 
where $ a and T a are the orthogonal fluxes 

$ a (r ; n, T) = (l — Pt) / a ( r . n> T; 0), T a (r ; n, T) = (1 - P) 7 a (r ; n, T), (130) 

with the adjoint projection operator P' 1 " given by Eq. (89). 

Hence, there is no constant component of the correlation function due to the invariants. 
Such a time independent part would not lead to a convergent limit for the time integral, 
as required for the transport coefficients. The property expressed by Eqs. (130) is similar 
to the presence of the "subtracted fluxes" in the Green-Kubo expressions for the transport 
coefficients of molecular fluids, with the subtracted fluxes being orthogonal to the global 
invariants of the dynamics. 
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VI. EULER ORDER PARAMETERS 



At Euler order, the unknown parameters of the phenomenological hydrodynamics of 
Section II are the cooling rate ( (n,T), the pressure p(n,T), and the transport coefficient 
( u associated with the expansion of the cooling rate to first order in the gradients. The 
cooling rate has been already defined by Eq. (50) above. The pressure and ( u can also be 
given explicit definitions in terms of the correlation functions from the coefficient of k in 
Eq. (110). Since ( u was previously identified as given given by Eq. (112), comparison of the 
remaining Euler coefficients in Eqs. (20) and (110) gives 



\ 



= k ■ D(n, T; 0, 0), 



(131) 



hyd,(a) 
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nm 
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(132) 



/ 



Using the definition in Eq. (98) and taking into account once again Eq. (65), it follows that 

-dp h (F;n,T, £/)" 



k ■ D a/3 (n, T; 0, 0) = V~ l J dTk- f a (T; n, T; 0) 



dyp 



(133) 



[7=0 



It is shown in Appendix E that all the matrix elements of (132) follow from (133) if the 
pressure is identified as 



p(n, T) = {Vd)- 1 J dTp h (r; n, T) tr H(r), 



(134) 



where tr H(r) is the volume integrated trace of the microscopic momentum flux. Its detailed 
form is given by Eq. (E5) of Appendix E. This is the second non-trivial result of the linear 
response analysis here, providing the analogue of the hydrostatic pressure for a granular fluid. 
It is possible to show that Eq. (134) leads to p(n,T) = nT in the low density limit, but 
at finite density the dependence on temperature and density of the pressure is determined 
by details of the HCS distribution, rather than the Gibbs distribution. This is in contrast 
to the appearance of Cq (n, T) in Eq. (131) which is a choice made in the definition of the 
temperature. In general, there is no relationship of p(n, T) to e (n, T) via thermodynamics, 
as for a normal fluid. 

The transport coefficient ( u represents dissipation due to inelastic collisions proportional 
to V • U. It has no analogue for normal fluids, where the Euler hydrodynamics is referred 
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to as "perfect fluid" equations, since there is no dissipation in that case. The simplest 
representation of ( u is the intermediate Helfand form, Eq. (118). More explicitly, it is 
shown in Appendix E that it can be expressed as 

C u (n, T) = WmiVTe^Td)- 1 J dTW (F; n, T) e'^M^u (F; n, T) , (135) 



with the source term W (F; n, T) defined by 

'd(e , T T( ) 



W(T;n,T) = —LE(F) — N 



dn 



~ E(T) 



The phase function A4^u (F;n,T) is the conjugate momentum 



M c u (F; n, T) = drr 



Sp, 



5U{r) 



N 



d(e , T T( ) 
de 



dp h (T;n,T) 



(136) 



(137) 



{j //3 }={n,T,0} s= i 

The second equality makes use of the local HCS distribution form for the velocity dependence 
at uniform density and temperature, 

[Plh (r| {yf3})] M ={n,T,U(r)} = Ph [{Qr, V r - U (q r )} J 71, T) . (138) 

The corresponding Green-Kubo form follows from a similar analysis of (112), or by direct 
integration of Eq. (135), 



C U (n,T) = UmiVTeo^d)' 1 J dTW (F;n,T) M c u (F;n,T) 



limiVTeo^d)' 1 / dt' / dF W (F; n, T) e' Crt ' CtM.u (T; n, T) . (139) 



This completes the identification of the Euler order parameters of the linearized hydrody- 
namic equations. Namely, exact expressions for the pressure and ( u in terms of correlation 
functions for the reference HCS have been derived. 



VII. NAVIER-STOKES TRANSPORT COEFFICIENTS 

The six transport coefficients at order k 2 can be easily identified in terms of elements of 
the correlation functions matrices A and Y introduced in Eq. (110). The twelve interme- 
diate Helfand and Green-Kubo forms are given in Appendix E. Only the shear viscosity is 
discussed here in some detail. Consider first its intermediate Helfand form. The analysis 
parallels closely that of ( u in Appendix E, with the result 

77 = - lim V^ 1 J dFH xy (F)e~ ZTt M v (F; n, T) . (140) 
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Here, Hy(r) is the volume integrated momentum flux of Eq. (E5) in Appendix F, and Mr, 
is the moment defined by 

N 



Mr 



,(r;n,T)= [drx\^-\ = — ^ q rx -^—ph(T; n, T). (141) 

The xy components occur here since the x axis has been taken along k and the y axis along 
the transverse direction e\, to simplify the notation. The corresponding Green-Kubo form 
is 

j) = -\imV~ 1 J dVH xy {V)M v {V]n,T) 

+ lim J dt'V- 1 J dTH xy (T)e- ZTt 'C T M v (T; n, T) . (142) 

In this case, the projection operators in Eq. (129) can be omitted since their contributions 
vanish from symmetry, and the dynamics is orthogonal to the invariants without such terms. 

The transport coefficient ( u vanishes for normal fluids, but the shear viscosity remains 
finite, as it is well known. It is instructive at this point to compare and contrast the results 
(140) and (142) for normal and granular fluids. Suppose from the outset a local equilibrium 
canonical ensemble corresponding to the equilibrium p c (r), had been used to generate the 
initial perturbations. Then, assuming non-singular conservative forces, 

N d N 
Mr, (r; n,T) = -J2 qr X -^Pc(T) = mT- l p c (T) ^ q rx v ry = T- l p c iY)M xy , (143) 

r=l rx r=l 

with 

N 

M xy = ^2q r ,xV r ,y (144) 
r=l 

Then, 

C T M V (r; n, T) = T~ l Pc {T)LM xy = T~ l p c {T)H xy (T). (145) 

Therefore, the intermediate Helfand and Green-Kubo expressions for a normal fluid then 
become 

77 = - limiVTr 1 (H xy M xy (-t)) c , (146) 

and 

77 = hm(\/T)- 1 / di! (H xy H xy (-t')) c , (147) 
Jo 
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respectively. The brackets denote an equilibrium canonical ensemble average, and the dy- 
namics is that of the Liouville operator L = L. These are the familiar results that have been 
studied and applied for more than forty years. 

To make the comparison between the elastic and inelastic cases, consider first the inter- 
mediate Helfand forms (140) and (146). The similarity between the structure for the normal 
and granular fluid results is striking, but the substantive changes are significant. For the 
granular fluid, the equilibrium ensemble has been replaced by the HCS ensemble. In addi- 
tion, the Liouville operator has been replaced by that including the nonconservative force 
and L ^ L. Finally, the generator for the dynamics includes the effect of temperature cool- 
ing in the reference HCS, L Ct = L — ( Td/dT. These differences manifest themselves 
in the Green-Kubo expressions in analogous ways. The inclusion of nonconservative forces 
implies a time independent contribution, the first term of Eq. (142), which vanishes in the 
elastic limit of Eq. (147). Also, due to the change in the ensemble, the two fluxes of the 
time correlation function differ for a granular fluid, while both are momentum fluxes for a 
normal fluid. Still the structure is such that these fluxes are orthogonal to the invariants of 
the dynamics in both cases, so that the time integrals can be expected to converge. 

VIII. DIMENSIONLESS FORMS AND SCALING LIMIT 

The analysis presented up to this point is quite general, and the only restriction placed on 
the nature of the microdynamics is that it be Markovian and the trajectories be invertible. 
These restrictions are satisfied by most models used to describe the interaction between 
granular particles. Examples are the Hertzian contact force model [32], the linear spring- 
dashpot model [32] , and the system of inelastic hard spheres with impact- velocity dependent 
coefficient of restitution [33]. The definition of these models and the associated generators 
of phase space dynamics are shortly reviewed in Appendix B. 

In general, there are two energy scales in the problem being addressed here. One is the 
total energy per particle or, equivalently, the cooling temperature T h (t). The other energy 
scale is determined by a property of the specific collision model, called e in the following. 
For the Hertzian spring case, it is the average compression energy of the spring. For hard 
spheres, it is fixed by some characteristic relative velocity in the dependence of the restitution 
coefficient on the relative velocity of the colliding pair. It is useful to reconsider the Liouville 
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equation in a dimensionless form that identifies these two different scales. In the limit that 
their ratio e/T h (t) is small, a special scaling form of the results above is obtained. This limit 
is exact for hard spheres with constant restitution coefficient (e = 0), and makes precise the 
conditions under which that idealized model may be approximately valid for many states of 
interest. 

Consider the Liouville equation (54), and introduce the dimensionless variables 

s = s(t,T), (149) 

where 

MT) . (-) (150) 

is a thermal velocity, I is the mean free path, and the function s(t,T) verifies the partial 
differential equation 

with the boundary condition s(0,T) = 0. The corresponding dimensionless distribution 
function is 

p*(r*;e*,s) = [lv (T)]) Nd p(T;n,T,t), T* = {q* r ,v* r }. (152) 

The dependence on the (dimensionless) density has been omitted on the left hand side to 
simplify the notation. In these variables, the dimensionless Liouville equation takes the form 

d sP * (r*; e*, s) + T (r* ; e*) p* (V; e* , s) = 0, (153) 

where 

— * dp* Cn(e*) d —* 

C* (r* ; e*) p* = Co (e*) e*^ + ^Li ^ _ . ( v ; p *) + £, (r ; e *) p*, (154) 

r=l r 

with I* (r*; e) = lL(T)/v (T) and C * = l( (T)/v (T). 

To interpret this result further, it is useful to consider the distribution of the HCS, 
p* h (T*; e*), which is the steady state solution of Eq. (153), i.e., 

r(r ;e *M(rv*) = o. (155) 

This solution is the dimensionless form of the universal function ph (T;n,T), where the 
dependence on T has been separated into a part that simply scales the velocities, and a 
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part that adimensionalizes the collisional energy scale e. This shows that in the appropriate 
variables, the distribution function of the HCS is stationary and universal, even when velocity 
scaling alone (see below) does not hold. For an isolated system, e*(t) = e/2T h {t) grows with 
increasing t since the system temperature decreases. For very large e*(t), the collisions 
become elastic and the system approaches a normal fluid. However, the alternative view of 
(153) is to specify the solution as a function of (r*,e*) and then study its properties as a 
special non-equilibrium steady state of granular fluids. 

As just noted, for large e* the collisions become practically elastic. In the opposite limit, 
e* << 1, the dependence on e* of the distribution function can be neglected, 

p*(r*;e*,s)^P*(r*;s) (156) 

and the Liouville equation becomes independent of e* 

d s p*(T*) + £*(r*)p*(r*) = 0, (157) 

r* N a 

c*(n P *(T*) = k £ JL . KV(r*)] + (r*). (iss) 

r=l r 

This is the limit in which all temperature dependence occurs through velocity scaling alone, 
as there is no other significant energy scale. It occurs for sufficiently hard interactions and/or 
sufficiently large kinetic energy. Many simplifications then occur. The HCS solution has the 
simple form p* h (T*) and consequently Co is a pure number. Equation (151) defining s can 
be integrated in this case to give 



4 in 

Co 



\ _ Qv (T)t 



(159) 



21 

Similarly, the dimensionless form of the cooling equation for T h (t), Eq. (47), can be integrated 
to get the explicit dependence on t 

r„(0) 1 + 21 I [ ' 



Finally, when (158) is evaluated at T h (t) the relationship of s to t becomes 

(161) 



Co 



i+ 21 Qot 



The conditions for which Eq. (156) applies will be called the "scaling limit". 

The special collisional model of inelastic hard spheres with constant restitution coefficient 
has no intrinsic collisional energy scale, so e* = and the scaling limit is exact. The 
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generators for this dynamics are indicated in Appendix B and can be understood as the 
singular limit of a soft, continuous potential, like the Hertzian contact force model. The 
analysis of the preceding sections is specialized to this case in the following paper [22], where 
it is shown that the above simplifications admit a more detailed exposition of the formal 
expressions for the transport coefficients. The rest of this section is a brief translation of 
some of the main results here to their dimensionless, scaling limit form. 
The dimensionless hydro dynamic fields are defined by 

<*>-{£}-{£•!■;&}• 

where the definition of the y a h 's follows from the second identity. Then, the fundamental 
linear response equation, Eq. (75), in the dimensionless variables <5y* is 

5y*(k*,s)=C*(k*,s)5y*(k*,0), (163) 

with 

C* ap {k\ s) = Lf.s. Cap K, T h (t); k, t) y^ h [T h (0)] . (164) 
It follows from Eq. (80) that C*^ (k*, s) obeys the equation 

[d s + K,*(k*,s)]C*(k*,s) = 0, C*(k*,0) = I, (165) 
the dimensionless transport matrix being 



vo[Th(t)]y a , h [T(t)]' 



Kl p (k*, s) = -6 a pp a CZ + - ]rr !ZL L^ Kap (n, T; k, t) . (166) 



Here /C (n, T; k, t) is the transport matrix analyzed in the previous sections. The additional 
contributions to /C*(fc*, s), proportional to {p a } = {0, 1, |, . . . , |} arise from differentiating 
the normalization constants with respect to T. Because of the scaling limit, JC* a/3 (k*,s) is 
independent of T(t) and the hydrodynamic limit can be identified as 

1C hvd {k*) =/C*(fc*,oo). (167) 

The phenomenological form of IC* hyd (k*), corresponding to that of Sec. II above, is given in 
the following companion paper [22]. 
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The dimensionless forms for the response functions C*^ (fc*, s) as phase space averages, 
follow from Eq. (77), 



J dTa a (T;n,T, k) e~ tCT ip^{T\ n, T, —k) 

dTa a [T; n h , T h (t); k] e^^T; n h , T h (0); -k)y^ h [T h (0)} 



Vy a , h (T) 
l 



yp, h [T(o)] 

n=n h ,T=T h (t) 



Vy a , h [T h (t)} 



= V*- 1 / dT* a* a (P; fe*) e- sC i/jUT*; -k*), (168) 



where V* = V/l d . More details of this transformation are given in Appendix B of ref. 
[22]. The generator for the dynamics C* is given by Eq. (157), and the dimensionless phase 
functions are 

a{ ' j i*y a>fc CO ' ( j 

^(r*;fe*) = K(T)f d ^(r;n,T;fc)^(T). (170) 
The hydrodynamic transport matrix }C* hyd (k*) is therefore given by 

/C*^(fc*) = /C*^ d (0) - lim | [9 S + /C*^ d (0)] C* (fe*, s)} C*- 1 (fe*, s) , (171) 

which is the dimensionless form of Eq. (81). In this way, all relation to the cooling of the 
reference state through a dependence on T has been removed. However, the dynamics of 
homogeneous perturbations of this state remains through 

C*(0,s) = e-^ hvd (°>. (172) 

The dimensionless correlation functions defining the transport coefficients are obtained 
in a similar way and have representations analogous to (168). An important difference is 
that the generator is C* — JC hyd *(0), indicating that the homogenous hydrodynamics is com- 
pensated. These simplifications and further interpretation are also deferred to the following 
paper. 

IX. DISCUSSION 

The objective of this work has been to translate the familiar methods of linear response 
for normal fluids to the related, but quite different case, of granular fluids. In both cases, the 
linear response to perturbations of a homogeneous reference state is described in terms of the 
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fundamental tools of non-equilibrium statistical mechanics. This microscopic formulation is 
then compared with the corresponding description from phenomenological hydrodynamics, 
and the unknown parameters of the latter are identified in terms of associated response 
functions. The analysis entails several steps, and at each stage there are technical and con- 
ceptual differences encountered for granular fluids that have been addressed in the preceding 
sections. 

The first difference is the form of the reference homogeneous state. For normal fluids, 
this is the equilibrium Gibbs state, while for granular fluids the corresponding role is played 
by the HCS. In both cases, the ideal global case considered is thought to be representative of 
more general cases, where these global states are expected to represent the state locally. The 
Gibbs state is strictly stationary and therefore constructed from the dynamical invariants. 
The HCS has an inherent time dependence due to collisional energy loss ("cooling"). Here, 
it has been given a stationary representation by including the granular temperature as a 
dynamical variable. In this form, the granular linear response problem becomes similar to 
that for a normal fluid, namely the response to the spatial perturbation of a homogeneous, 
stationary state. However, the generators for that dynamics are now more complex, due 
both to the non- conservative forces responsible for collisional energy loss and a generator 
for changes in the granular temperature. 

The response functions contain information about both hydrodynamics and microscopic 
collective excitations. Therefore, the initial perturbation will excite in general a wide range 
of dynamics and the extraction of the hydrodynamic branch at long wavelengths, long times, 
can be quite complex. However, if the perturbation is chosen to excite only the hydrodynamic 
modes, this analysis becomes simpler and more direct. Practically, this can be done only 
in the long wavelength limit. For a normal fluid, these are the dynamical invariants, since 
the hydrodynamic modes are those that vanish in that limit. The associated perturbation 
is given by the corresponding global conserved quantities (number, energy, momentum). 
The granular fluid is more complex, since there is a residual hydrodynamics even in the 
long wavelength limit. This is due to the nonlinear temperature cooling that is linearized 
about a reference cooling state. The hydrodynamic modes are therefore identified from this 
non-zero long wavelength dynamics. The microscopic perturbations corresponding to this 
dynamics are described in Sec. III. Furthermore, by extracting this dynamics from the the 
microscopic evolution, these perturbations become the invariants of the residual dynamics 
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(see Eqs. (62) and (63)). 

With this knowledge of the special long wavelength hydrodynamic perturbations, spatial 
perturbations are constructed from their corresponding local forms. In the normal fluid 
case, these are the microscopic local densities of number, energy, and momentum, and 
are generated by a local equilibrium ensemble. By analogy, the densities of the invariants 
for the granular fluid are the appropriate perturbations, and they are generated from a 
corresponding local HCS. Since the HCS is not the Gibbs state, these densities are no longer 
the local conserved densities for a normal fluid. 

The remaining analysis is straightforward, extracting the Euler and Navier-Stokes pa- 
rameters as coefficients in a wavevector expansion, using the conservation laws to assist in 
the ordering. There are two different sets of conservation laws (balance equations) for the 
granular fluid, in contrast to only one set for the normal fluid. Consequently, expressions for 
transport coefficients are not simply given in terms of autocorrelation functions of appropri- 
ate fluxes. Instead, for the granular fluid there are two conjugate fluxes in the Green- Kubo 
expressions. 

The utility of these results rests on further studies of appropriate ways to evaluate them. 
The status now is that exact expressions for the quantities of interest (e.g., pressure and 
transport coefficients) are given formally without any inherent uncontrolled approximations 
(e.g., as in some chosen kinetic theory). This is a more suitable point for the introduction 
of practical methods for evaluation. In the study of normal fluids, a number of methods 
have proved very useful. They include molecular dynamics simulations, memory function 
models incorporating exact initial dynamics, and linear kinetic theory. It is hoped that 
similar approaches will be developed for the expressions provided here. To elaborate on 
this, consider again the new Euler order transport coefficient given by Eq. (135), make the 
temperature scaling of the generator explicit, and evaluate the entire expression at T = T h {t) 



t u [n,T h (t)] = -lim 



(VTeo^d)- 1 [dTW (T; n, T) e^V^ jh q s ■ ^i^T) 



s=l 

N 



= -Hmtyr^fKTpThWlC 1 / dTW[T;n,T h (t)]e- Lt J2q s - 

J s=l 



T=T h (t) 

dp[T;n,T(0)] 
dv s 

(173) 

The generator for the dynamics is now that for the trajectories alone, so this is a form 
suitable for MD simulation. The simulation method must be constructed in such a way as 
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to represent the unknown HCS appearing in (173). Elsewhere, the evaluation by kinetic 
theory [34] is considered, and the original representation given by Eq. (135) is found to be 
more suitable. 

It is worth recalling that liquid state transport for simple atomic fluids remains a pro- 
totypical strongly coupled many-body problem, with little progress beyond simulation of 
formal expressions such as those given here. More should not be expected for "complex" 
granular fluids. The formal representations of transport coefficients by methods of statisti- 
cal mechanics provides a new perspective on a difficult old problem. As for normal fluids, 
significant further progress can be expected for the idealized model of hard spheres. That 
is the subject of the following companion paper. 
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APPENDIX A: HOMOGENEOUS STATE DYNAMICS 

The dynamics associated with the homogeneous cooling state of interest here is two fold. 
The first is the cooling of the temperature, determined from the solution to Eq. (15). For a 
given initial condition T, the solution is denoted by T h (t;n h ,T). The density is a constant 
parameter and sometimes it is left implicit in the notation of the text. The second dynamics 
is the linear response to small homogeneous changes in the initial conditions, 

( dT h (t-n h ,T) \ f dT h (t;n h ,T) \ 
6T h (t, n h ,T)= ^ — jjn h+ ^ j ^ 6T. ( Al) 

In this Appendix, it is shown how the response function for this second type of dynamics is 
obtained from the linearized hydrodynamic equations to give Eq. (33). 
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A useful identity for any function of the temperature, X(T) is 

X[T(t 2 )} = exp |- (t 2 - tO Co [T (h)] T ( tl ) } X[T( tl )]. (A2) 

Here T(t) is a solution to Eq. (15). The identity can be proved by performing a Taylor 
series of X[T(t 2 )] in powers of (t 2 — ti) and using Eq. (15) to evaluate the time derivatives 
in terms of T(ti) derivatives. The above identity gives, in particular, 

d 



T(0) = e W ^tCo[T(t)]T(t) dT{f) 
T{t) = exp(-tCo[T(0)]T(0) 



and 



( dT(t) \ 



dT(0) 
(o[T(t)]T(t) 



T(t), 
T(0), 



(A3) 



(A4) 



\dT(0)J nh Co[T(0)]T(0) 
Consider some function X[T(t;T),t] that depends on time through T(t;T) plus some 
residual time dependence. Use the second equation of (A3) to write 

d 



X [T(t;T),t] = exp 
and, consequently, 

d t X[T(t;T),t] = exp -t( (T) T 

d t -( (T)T 



-Ko(T)T 



dT 



X(T,t) 



(A5) 



d_ 
dT 
d_ 
df 



d t -( (T)T 



d_ 
df 



X(T,t) 



X(T,t) 



(A6) 



T=T(t;T) 

The time dependence due to T(t; T) can be replaced by treating T as an independent 
variable, with the additional generator for its dynamics Co (T) Td/dT. It is then equivalent 
to determine X(T,t), and evaluate it finally at T = T(t;T). In the case of Eq. (26), this 
leads to 

C hyd (n,T;k,t) = 0, (A7) 



dt - Co (T) T— + K, hyd (n, T; k) 

Caf (n,T;k,0) = 5 a p, 
with the definition in Eq. (29). For k = 0, use of Eqs. (20)-(22) gives 

^(n,T;0,t) = ^, 



for a 7^ 2, while for a = 2 it is found: 



dt-Co(T)T-^+]C^ d (n,T;0) 



C\f (n, T; 0, t) + K,T (n, T; 0) S w = 0, 



(A8) 



(A9) 



(A10) 
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or, more explicitly, 



d 

d t - CoT— + 



d (CoT) 
dT 



C*f(n,T;0,t) + 



The solution of this equation is 
C%f(n,T;0,t) 



dn J 



5 



T(-t;T) 



1/3 



d(CoT) 
dn 



dT 



'1/3 



(All) 



dT(-t;T)J n 



^2/3, 



(A12) 



as can be verified by direct substitution into Eq. (All) and repeated use of Eq. (A4). This 
is the result (33) of the text. 



APPENDIX B: GENERATORS OF DYNAMICS 



The interaction between the constituent particles of the dissipative fluid enters the pre- 
sentation here via the Liouville operators that generate the dynamics. The analysis of the 
text places few restrictions on these generators and admits a large class of models to repre- 
sent real systems. For example, it is not necessary that they be pairwise additive, although 
the examples of this appendix all assume that case. There is a qualitative difference between 
the generators for continuous or piecewise continuous forces, and those for singular forces 
(e.g., hard spheres). Examples of each are given here for illustration. 



1. Dissipative soft spheres 

The fluid is assumed to be comprised of mono-disperse spherical particles with pairwise 
additive central interactions. The latter implies that the forces are "smooth", without 
tangential momentum transfer, and Newton's third law holds. The simplest realistic model 
for the force F that particle s exerts on particle r is the smooth, frictional contact model 
[32, 33] given by 

F (q ra , g rs ) = q rs <d (cr - q rs ) [f (a - q rs ) - 7 (a - q rs ) (g rs ■ q rs )] . (Bl) 

where O(x) is the Heaviside step function, q rs = q r — q s is the relative coordinate, g rs = 
v r — v s is the relative velocity of the two particles, and q rs = q rs /q r s is the unit normal 
vector joining the centers of the two particles. Moreover, f(x) and 7 are a function and 
a constant, respectively, to be described below. This is a piecewise continuous force that 
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vanishes for separations greater than a, which therefore can be thought of as the diameter of 
the particles. The first term between the brackets describes a conservative force representing 
the elastic repulsion due to the deformation of real granular particles. If / (x) is chosen to 
be linear, the deformation is that of a spring. The amount of deformation can be adjusted 
by the choice of the spring constant. A second, more realistic, choice is the Hertzian contact 
model for which / (x) oc x 3 ^ 2 . 

The second term of (Bl) is a nonconservative force representing the energy loss of the 
particle pair on collision. It is proportional to the relative velocity of approach during the 
collision, and the amount of energy loss is adjusted by the choice of the friction constant 7. 

The Liouville operators L and L, defined in Eqs. (37) and (38), for the dynamics of phase 
functions and distributions for these models can be identified as 

N 8 1 N N 8 

LX (r)^J2 v r .g-X(T) + -J2J2 F («™, dr.) ■ ^ (F) (B2) 

r=l yr r=l s^r r 

and 

1 N N 8 

LX (r) = LX (r) H ^2 E X ^7T- ■ F (I™, Qrs) ■ (B3) 

r=l rj±s 

It is readily verified that the total momentum is conserved, since Newton's third law is 
satisfied, i.e., F (q rs ,g r s) — —F (q S r,9sr)- The total energy is 

AT N N 

e(t) = + - E v (*■•) ' ( B4 ) 

r=l r=l s^r 

where the potential energy function V(^) verifies 

= _ e (a - q rs ) f(a- q rs ) . (B5) 

dq rs 

The microscopic energy loss is easily computed by using Eq. (40) with the result 

^ N N 

LE(Y) = -- ^ E 6 (°" ~ qrs "> 7 (°" _ ( 5rs ' ' ( B6 ) 

r=l s^r 

showing that it is associated with the nonconservative part of the interactions as it should. 



2. Hard sphere dynamics 

For a given energy of activation, the contact forces considered above may have small 
deformations, i.e. the region in which the forces differs from zero verifies (a — q rs ) /er <C 1. 

44 



In that case, the conservative part of the force approaches that of elastic hard spheres. The 
primary effect of the nonconservative force is to decrease the magnitude of g rs ■ q rs after the 
collision. This can be represented by the scattering law 

9'rs = 9rs ~ [1 + « (9rs)] " 9rs) (B7) 

where g' rs is the relative velocity after collision and a (g rs ) is a coefficient of restitution 
that depends on the relative velocity. The total momentum of the pair is, by definition, 
unchanged in the collision. The elastic limit corresponds to a (g rs ) — > 1. Subsequent to 
the change in relative velocity for the pair (r, s) , the free streaming of all particles continues 
until another pair is at contact, and the corresponding instantaneous change in their relative 
velocities is performed. The collision rule is assumed to be invertible, i.e., a (g rs ) is specified 
so that the trajectory can be reversed. 

Since there is no longer a potential energy, the total energy for the system is its kinetic 
energy, which changes on a pair collision by 



)m (g' r 2 s - g 2 rs ) = -\m [l - a 2 {g rs )] (a ■ g rs f . (B8) 



4 v rs ^ rs ' 4 

This is clearly the analogue of the terms on the right hand side of Eq. (B6). In fact, the 
velocity dependence of the coefficient of restitution can be modeled from a comparison of 
the two equations. 

There are two components of the generators L and L, corresponding to each of the two 
steps of free streaming and velocity changes at contact, 

N o 1 N N 

L = H v r-^r + -Y,Y, T ^sl (B9) 

r=l r r=l s^r 

N „ . JV JV 

T =H v r-^r-- 2 Y.Hnr,s). (BIO) 

r=l r r=l s^r 

The operators T(r, s) and T(r, s) describe the binary collision for a pair, 

T(r, s) = 5(q rs - a)Q(-g rs ■ q rs )\g rs ■ q r s\(b rs - 1), (Bll) 

T(r, s) = 5(q rs - a) [j (v r , v s ) b^ 1 - l] <d(-g rs ■ q rs )\9rs • q ra \- (B12) 
Here b rs is a substitution operator, 

b rs X(g rs ) = X(b rs g rs ) = X(g' rs ), (B13) 
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which changes the relative velocity g rs into its scattered value g' rs , and b r } is its inverse. 
Finally, J (v r ,v s ) is the Jacobian for the transformation from {v r ,v s } to {v' r ,v' s }, 

J (v r ,V a ) = 

The delta function in (Bll) and (B12) requires that the pair is at contact, while the theta 
function requires that the directions of velocities are such as to assure a collision. A deriva- 
tion of these results and further details are given in the companion paper following this 
one. 



d (b rs v r ,b rs v s 



d(v r ,v s ) 



(B14) 



APPENDIX C: HOMOGENEOUS COOLING SOLUTION 

In this Appendix, Eqs. (58) and (59) are proved, leading to the exact solution of the 
Liouville equation (60) for homogeneous perturbations of the HCS. The HCS is the stationary 
solution to the Liouville equation (54), i.e., 

CTp h (T;n,T,U)=0, (CI) 

Z T = -( (n,T)T-^ + L. (C2) 
The action of the operator C T on \P a can be evaluated as follows: 

= ^^(r;n,T,J7)/C^(n,T;0), (C3) 

P 

where /C^f(n, T; 0) has been identified in the last line from Eqs. (20)-(22) particularized for 
k = 0. This proves Eq. (58). 
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Next, using this result, 
C 2 V a (T; n, T,U) = J2 pr*/^; n, T, U)} /Cg d (n, T; 0) 
+ */»(r; n, T, U)£ T KX(n, T; 0) 
= E E ^(r; n, T, C7)/Cy (n, T; 0)/C^(n, T; 0) 

(3 7 



+ J]* /3 (r;n,T, C7) 

= EE^(r;n,T,c/) 

/3 7 



<o(n,T)TA 



/C^(n,T;0) 



(n,T; 0) - * 7/J Co(n, T)T 



<9T 



/C?>,T;0), 
(C4) 



or, in a compact matrix notation, 



L 2 T ^ = V ( K, hyd - KoT-^j JC hyd , 



(C5) 



with 7 denoting here the d + 2 unit matrix. Applying the induction method, 



= (CtV) ( /C^ - /CoT^ ) JC hvd + ^£ 



<9T 



1-2 



fchyd _ 



dTj 



^K hyd ) (/C^ - ^CoT^ /C^ + * (-CoT^) (/C fc * - /CoT^^) ^ 



= * ^/C^ - KvT—j K, hyd 
This implies 

e- Z ^ a (r;n,T,*7) = E*/ 3 G"; n, T, E7) jexp 



(C6) 



^(n, T; 0) - /Co (n, T) T A 



= ^ * (r; n, T, C7) Cjf (n, T; 0, t) , 



/3a 

(C7) 



which proves Eq. (59). In the last transformation, the formal solution of Eq. (30) has been 
used. 
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APPENDIX D: MICROSCOPIC CONSERVATION LAWS (BALANCE EQUA- 
TIONS) 



1. Fluxes associated with a a (T;n,T;k) 



The microscopic balance equations for the phase functions a a (T; n, T; k) follow from those 
for the Fourier transformed number density J\f(T; k), energy density £(T; k), and momentum 
density G(T; k) defined in Eq. (68). These balance equations relate the time dependence of 
the densities to appropriate fluxes 

/ g(r ; fc) \ / n \ 



d t e 



Lt 



( N{T-k)\ 
£{Y;k) 



— ik • e 



Lt 



s(T;k) 

\h(r-k) ) 



-Lt 





w(T;k) 




(Dl) 



These are microscopic conservation laws for Af (T; k) and Q (r; k). For granular fluids, the 
energy density has a source w (T; k) due to the inelasticity of the collisions. The forms of 
the fluxes of M (T; k) and Q (k) are obtained from 



LM (r ; k) = — k ■ Q (r ; k) , LQ (r ; k)=ik-h (T; k) . 

m 



(D2) 



The expression for the tensor momentum flux h is 

N -i N N 



hii (r; k) = mvr^rje*-* + - dxJ^Yl 0™) e ik <»*° + *\ (D3) 

r=l ^° r=l s^r 

This is the usual result for nonsingular forces F, generalized here to include a noncon- 
servative contribution as well. Some examples are discussed in Appendix B. In all of this 
Appendix only nonsingular forces are considered. The corresponding results for hard spheres 
are given in the following companion paper. 

The right sides of Eqs. (D2) are proportional to k, indicating that they are densities of 
conserved variables. For the energy density there is both a flux and a source 



L£ (T; k)=ik- s(T; k) - w (T; k) . 



(D4) 



The energy flux is given by 



N 



s(T;k) = 



r=l 



mv 
~~2 



2 i N 



s^r 
N N 



ikq r 



+ - / dx qrs ( Vr + V ^' F ^rs, 9rs) e tk 

^° r=l s^r 



(xq rs +q s ) 



(D5) 
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and the source term is 

1 N N 

w(T;k) = --J2Y,9rs- F nc (q rs , g rs ) e lk (D6) 

r=l s^r 

where F nc (q rs ,g rs ) is the nonconservative part of the force. The functional forms for the 
fluxes G (T; k) , s(T; k) , and h (I"; k) are the same as those for a normal fluid, except that 
the total force, including its nonconservative part, occurs. The source w (T; k) depends only 
on the nonconservative part of the force. For the special case of the force given in Eq. (Bl), 
Eq. (D6) becomes 

N N 

w(T;k) = -J2J2 e ^- 9") 7 ( a " 9") ' «") 2 e lk q \ (D7) 

r=l s^r 

which agrees with Eq. (B6) for = 0. 

The corresponding fluxes associated to the a a (T; n, T; k) follow from their definition, Eq. 
(67), in terms of the above densities, 

La a (T; n, T; k) = ik ■ f a (T; n, T; k) — w(T; k)5 a2 , (D8) 

eo,T 

with _ 

/ 1 (r;n,T;fc) = ^ff^, (D9) 
m 

f 2 (T; n, T;k) = — \s(T; k) - (r ; fc)l , (D10) 

e ,T L rn J 

\ l3 (T-,n,T-k) = ^^. (Dll) 

The last equation above gives the tensor flux associated to the vector a 3 = Q /nm. 

Next, calculate the quantity £(T;n,T; k) appearing in Eq. (85). Use of Eq. (84) gives 
directly Eq. (85) with 

I(T; n, T; k) = —w (T; k) + Co (n, T) T JU 2 (r ; n, T; k) - V /Cg d (n, T; 0) a^r; n, T; fc) 

1 ^ Co(n, T)T ,9eo,T~ rr T Co(n,T)T de , n „ 
-w (T; k) ^-^a 2 (T; n, T; k) ^-a^F] n, T; k) 



e ,T e ,T dT e ,T dT 

-ai(T; n, T; k) — — -a 2 {T; n, T; k) 



On 1V ' ' ' ' dT 

1 f 3{T;k) _( °Mn,T)n ) , l(r; „ iT;fe) 



eo,T 1 \ <9n 

'9[e ,TCo(n,T)T] 



dT 
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a 2 (r>,T;fc) \. (D12) 



From the expression for the cooling rate in the HCS given in Eq. (50), 

e ,TCo(n, T)T = -V' 1 J dTp h (T; n, T)LE{T) = V' 1 J dTp h (T; n, T)w (T; 0) (D13) 
and, since w (T; 0) is independent of n and T as seen from its definition in Eq. (Dl), 

{ d{e °dn T) ) T = rfr ^( r ^' T )-( r '°)' 

( d{e ° d f T) ) =V- 1 /dT* 2 (r;n,T)^(r;0). (D14) 

Substitution of the above relations into Eq. (D12), noting that the sum can be extended 
to include a — 3, since the new contribution vanishes by symmetry as a consequence of w 
being a scalar, gives Eq. (86) in the text. 

2. Fluxes associated with ip a (F;n,T;k) 

By definition in Eq. (65), the set of ip a (T;n,T;r) are densities associated with the in- 
variants, i.e. 

^(r>,T;0) = ^ a (r>,T), (D15) 

and so 

U a p (t, T) MT; n, T; 0) = MT; n, T; 0) (D16) 

P 

and 

d t MT;n,T;0,t)=0. (D17) 

The time derivative of ip a (T;n,T; k,t) must be of order k, so there exists a flux 
7/3 (r; n, T; k, t) such that 

dtMT; n, T; k, t) - ik-% (r ; n, T; k, t) = 0. (D18) 

The generator of the dynamics U (t, T) is defined by Eq. (63), and taking the time derivative 
there, it is seen to obey the equation 

d t U (t, T) = — [C T - IC hydT (n, T; 0)] U (t, T) , (D19) 

where K, hyd T is the transpose of JC hyd . This equation can be formally integrated to write 

U(t,T) = exp{-t \£ T -lC hydT (n,T;0)]} , (D20) 
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that shows that Eq. (D19) is equivalent to 



d t U (t, T) = -U (t, T) \Z T - JC hyd T (n, T; 0)] 



(D21) 



Use of this into Eq. (95) yields 



d t %l)(T- n, T; k, t) = -U{t, T) \C T - JC hyd T (n, T; 0)] i/j(T; n, T; k). 



(D22) 



Comparison with Eq. (D18) leads to the identifications given in Eqs. (101) and (102). 



APPENDIX E: DETAILS OF EULER ORDER PARAMETERS 



The Euler order time independent correlation function k-D a p(n, T; 0, 0) whose expression 
is given in Eq. (133), is determined from direct evaluation. Consider first the case a — 1 for 
which /i(r ; n, T; 0) = Q (r ; 0) /m. Then 

-dp h (T;n,T, U) 



k- D lf5 (n,T;0,0) = {mVy 1 / dVk ■ Q{T; 0) 

d 



[7=0 



{mvy 1 



dTk-g(T;0) Ph (T;n,T, U) 



u=o 



v- 1 



d 



Nk ■ U 



u=o 



(El) 



in agreement with Eq. (131). For a — 2, use Eq. (D10) to get 
k-D 2(3 (n,T;O,0) = (Ve^y 1 f dVk ■ \s{V- 0) - ^g(F; 0) 



dp h (T;n,t,U) 



u=o 



d 



-e , n (mV 'eo.r) -1 



drk-s{T;0)p h {T;n,T, U) 
' d 



dyp 



u=o 

dTk-g(T;0) Ph (T;n,T, U) 



(E2) 



[7=0 



The second term on the right hand side is easily evaluated using the result in Eq. (El). 
The ensemble average in the first term can be carried out by making the change of velocity 
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variables v r — > v r + U and using the properties of the Galilean transformation, 
J dT s(T; 0)p h (T; n, T,U) = J dT [s(T; 0)} {Vr ^ Vr+u} p h (T; n,T,U = 0) 



dT { s(T; 0) + 



1 



S(T;0) + \mNU 2 



c/ + h(r ; o) • u 



+ -Q(T; 0)U 2 + UU ■ g(T; 0) J p h (r ; n, T) 
e (n,T)\/C/ + imiVf/ 2 C/ + J dTh(T;0) ■ U p h (T;n,T). 



(E3) 



Using this result it is easily obtained 

5/33 



k- D 2/3 (n,T;0,0) 



e ,r 



e - e^ n n + (Vrf)- 1 ^ dTtr H(T)p h (T; n, T) 



(E4) 



Here H(T) = h(r,0), so tr H(r) = Ylt=i 0) is the volume integrated momentum flux. 
For the case of the dissipative hard spheres discussed in Appendix B, it follows from Eq. 
(D3) that 



N ^ N N 

H *J ( r ) = ^2 mV r,i V r,j + 9 X] (<?r S , ■ 



(E5) 



r=l 



For a = 3, use Eq. (Dll) to get 

k-D 3p {n,T; 0,0) = (nm7)' 



g 

dy/3 



drh|m(r;0)p h (r ; n,T, t/) 



c/=o 



= (nmV)- l i^-j dT [h||||(r;0) + 2C/||^||(r;0) 
+mnU 2 ] p h (T;n,T,U = 0)} u=Q , 



(E6) 



where again the change of velocity variables has been made. The subindex || indicates the 
component in the direction of k. Therefore, 



k-D 3f) (n,T- 0,0) = (mn)- 1 (s^ + S^-^j (Vd)' 1 j dTti W(T)p h (T; n, T) . (E7) 

These results are consistent with the form of the phenomeno logical matrix ](j^ yd ^ i n Eq. 
(132), if the pressure is identified as 



p(n, T) = {Vd)' 1 J dT trH (r) p h (T; n, T) . 



(E8) 



This is Eq. (134) in the main text. 
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The single transport coefficient at Euler order ( u (n,T) was identified in Eq. (118) as 

T( u (n,T) = -limfe-5g(n,T;0, (E9) 
where the expression for S$ follows from Eq. (99), 

-V- 1 J dTl(T;n,T;0)^\T;n,T;t). (E10) 
The first term on the right hand side vanishes since 



/ 



dr£^(T;n,T)^ 3 (T;n,T) 



d 



dT£^(T;n,T)p h (T;n,T, U) 



(Ell) 



[7=0 



The same change of variables v r — > v r + U as above, shows that this vanishes due to spherical 
symmetry of p h (T;n,T, U = 0). The remaining contribution to ( u (n,T) as given by Eq. 
(E9), is made more explicit using the definition of £(T;n,T; 0), Eq. (88), and also that of 
U 3a (t,T), Eq. (63), to compute 

fii\T;n,T;t) = £ U 3a (t, T)^\r ; n, T). 



(E12) 



This gives 

( U (n,T) = +lim(VT)- 1 J2 [ drT(T;n,T;0)U 3a (t,T)k ■ ^(T;n,T) 

= -lim^Teo^)- 1 ^ f dT [(l - pt) L £(r)] S^e'^k ■ ^(F; n, T) 

= WmiVTe^Td)- 1 J dTW (T; n, T) e~ ZTt M c u (F; n, T) . (E13) 

The phase function W (T; n, T) is defined by 

W(T;n,T) = - (l - P f ) LE(T) 

= —LE(T) ~n(J^ (e ,TTCo)) - E(T) (eo^Co)) ■ (E14) 

The second equality follows from explicitly evaluating the action of the projection operator 
P', defined in Eq. (89), on LE(T). Finally, the phase function Ai(u (T;n,T) is 

M c u(T;n,T) = dk • (T;n,T)d 



— d J drr\\ 


5 pm 




m\ (r). 




N 
r=l 


dp h (T;n 


,T) 
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(E15) 



In the last transformation, the local equilibrium form for the velocity dependence has been 
taken into account, 



[Plh (r| {yp})]{y }={n,T,U{r)} = Ph [{Qv, V r ~ U (<? r )} ] U, T] 



(E16) 



APPENDIX F: NAVIER-STOKES ORDER TRANSPORT COEFFICIENTS 

The Helfand forms for the transport coefficients are identified from (119) and (120). For 
the energy flux, these are the thermal conductivity A and the new granular fluid coefficient 



u = eo,T lim kk 



d£ (n, T;t)-J2 D2a(n, T; 0, 0)Cg (n, T; t) 



D^(n, T;t)-J2 D 2a (n, T; 0, 0)C^(n, T; t) 



7^(1), 



(Fl) 
(F2) 



The shear and bulk viscosities, 77 and k, are identified as 



K + 



i] = mn lim kk : 
2(d- 1)77 



D^(n, T; i) - £ £> 4a (n, T; 0, 0)c£(n, T; f) 



(F3) 



mn lim kk 



(n, T; t) - ]T r> 3 a(n, T; 0, 0)C™ (n, T; i) 



(F4) 



Finally, the two Navier-Stokes transport coefficients associated with the cooling rate are 



C n = T~ 1 lim kk : S% (n,T;t) + T( u (n,T)k ■ C£ (n,T;t) 

( T = T" 1 lim [kk : Sg } (n, T; t) + T( u (n, T)k ■ cff (n, T; f) 
The corresponding Green-Kubo forms are 



A = e 0tT kk : 



D^(n, T; 0) - lim / dt' G 22 (n, T; t' 



(n, T; 0) - lim / dt' G 21 (n, T; 



K + 



n = nmkk : D$(n, T; 0) - lim / dt' G 44 (n, T; t') 

Jo 

2{ ~ ^ = nmkk : log (n, T; 0) - lim f dt' G 33 (n, T; t') 

Jo 



d 



(F5) 
(F6) 

(F7) 

(F8) 
(F9) 
(F10) 
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c = 



c 1 = 



T^kk : jsg^TjO) 

- lim J dt' [N^ (n, T; t') + T( u {n, T)kE 31 (n, T; 0, t') 
T- l kk-.[s { £{n,T-Q) 

-lim^ dt' ]^{n,T;f) + T( u {n,T)kE 32 (n,T;0,t' 



(Fll) 



(F12) 
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